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Foreword 


The papers presented here are from the Langley Research Center Workshop on 
Computational Fluid Dynamics (CFD) Validation of Synthetic Jets and Turbulent 
Separation Control (nicknamed “CFDVAL2004”), held March 2004 in Williamsburg, 
Virginia. The goal of the workshop was to bring together an international group of CFD 
practitioners to assess the current capabilities of different classes of turbulent flow 
solution methodologies to predict flow fields induced by synthetic jets and separation 
control geometries. This workshop consisted of three flow-control test cases of varying 
complexity, and participants could contribute to any number of the cases. Along with 
their workshop submissions, each participant included a short write-up describing their 
method for computing the particular case(s). These write-ups are presented as received 
from the authors with no editing. Descriptions of each of the test cases and experiments 
are also included. 
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CASE 1: Synthetic Jets in Quiescent Air 

C. S. Yao, F. J. Chen, D. Neuhart, and J. Harris 


Flow Physics & Control Branch , NASA Langley Research Center , Hampton , VA 23681-2199 


Introduction 

An oscillatory jet with zero net mass flow is generated by a cavity -pumping actuator. Among the 
three test cases selected for the Langley CFD validation workshop to assess the current CFD 
capabilities to predict unsteady flow fields, this basic oscillating jet flow field is the least complex 
and is selected as the first test case. Increasing in complexity, two more cases studied include jet 
in cross flow boundary layer and unsteady flow induced by suction and oscillatory blowing with 
separation control geometries. 

In this experiment, velocity measurements from three different techniques, hot-wire 
anemometry, Laser Doppler Velocimetry (LDV) and Particle Image Velocimetry (PIV), 
documented the synthetic jet flow field. To provide boundary conditions for computations, the 
experiment also monitored the actuator operating parameters including diaphragm displacement, 
internal cavity pressure, and internal cavity temperature. 

Experiment Setup 
Synthetic jet actuator 

The synthetic jet actuator for this experiment was based on an earlier design studied in detail by 
Chen et al (2000) 1 and Chen (2002) 2 . The jet exit slot of the actuator (Fig. 1) is 0.05 inches wide 
and 1.4 inches long. The slot width was larger than the previous design to allow velocity profiles 
across the slot to be properly resolved. The actuator is flush mounted on an aluminum plate at the 
floor, and is covered by a glass enclosure, 2x2x2 feet in dimension and l A inch thick. The glass 
enclosure isolates the synthetic jet flow from the ambient air and serves to contain the seeding 
particles as well. A sliding door on the glass panel provides access for seeding and the hot-wire 
probe. The jet is located at the center of the floor, and the slot is parallel to the glass sidewall. 

The actuator has a Plexiglas housing with a narrow cavity beneath the slot. The jet flow was 
driven by a single piezo-electric diaphragm, 2" in diameter, mounted on one side of a narrow 
cavity. The piezo diaphragm was driven at 445 Hz, which was selected to operate away from the 
cavity resonant frequency near 500 Hz. Figure 2 shows the output of the actuator as a function of 
the forcing frequency. An o-ring seal, 1.85" in diameter, is clamped between the diaphragm and 
the actuator cavity. Maximum jet velocity generated could reach 30 ~ 50 m/s, depending on the 
actuator performance, and the aging of actuator. In daily test, there was no performance 
degrading with time. The actuator did show a gradual decreasing of the pressure level with a long 
period of usage over months. Sometime, the diaphragm needed to be replaced due to the cabling 
faults. 

Actuator parameters 

Three actuator operating parameters including diaphragm displacement, cavity pressure and 
cavity temperature, were measured to monitor the actuator performance, and provide boundary 
conditions for computation and actuator modeling (Fig. 3). A fiber optics displacement sensor is 
installed to measure the diaphragm displacement at the center of the piezo-electric disc. The fiber 
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probe was calibrated in-situ using a micrometer. The displacement measured ranges between 
+0.2 ~ -0.4 mm, for in-ward and out-ward displacements. A dynamic pressure transducer 
installed at the center of the sidewall (the fixed wall) opposite of the diaphragm to measure the 
instantaneous cavity pressure waveforms. A thermal couple device is installed at the bottom of 
the cavity to monitor the internal temperature. 

Hot-wire measurements 

Hot-wire Anemometry provides a single point in space and a time history measurements of the 
velocity component perpendicular to the sensing wire. In the current setup, hot-wire 
measurements were made with a constant-temperature anemometer (CTA) and a single-sensor 
hot-wire probe. The sensing element of the hot-wire probe, 5 jam in diameter and 1.25 mm long, 
is a platinum-plated tungsten wire operated at an overheat ratio of 1.8. The probe was mounted 
perpendicular to the floor with the sensing element located at the center of the slot and parallel to 
the long axis of the slot. Measurements were taken at 47 stations with height ranges to 50 mm 
above the slot. A hot-wire probe support was mounted on a translation stage with a computer- 
controlled traversing system installed on the exterior of the glass enclosure. A slotted opening on 
the sidewall provides the access into the glass enclosure. 

The hot-wire probe was calibrated with a commercial desktop calibration unit. The signals 
from the CTA outputs were sampled at 100 kHz with 16384 data points recorded at each station. 
For a forcing frequency of 445 Hz, there are 72 periods of waveform recorded. Within each 
period of the driving cycle there were 225 data points. Seventy-two samples were averaged to 
compute the phase-locked statistics. Hot-wire signals were de-rectified to obtain negative 
velocities during the actuator suction cycle at lower stations. Uncertainty was +/- 2.5 % for 
measurements of heights of 5 mm or higher and +/- 6 % for lower stations near the solid surface. 

Laser Doppler Velocimetry measurements 

Similar to hot-wire measurements, LDV is single point in space measurement technique. LDV 
signals provide a velocity record and the corresponding particle arrival time whenever a valid 
signal is detected. In this test, a LDV system was set up to measure the vertical jet velocity and 
its decay along the centerline of the slot. Vertical velocity measurements were made at 47 
locations above the rectangular slot between 0 to 70 mm. The transmitted laser beams from a 
fiber optic probe were projected through the glass sidewall in a direction parallel to the long axis 
of the rectangular jet slot, and in a vertical plane. The crossing point of the two laser beams, 
which is the LDV sample volume, was centered over the center of the slot, both longitudinally 
and laterally. The scattered light from the seeding particles was collected by a second fiber optic 
probe, which was located in the forward-scatter direction, 50 degrees off of the direct forward 
scatter axis. 

The seeding particles were 0.9-micron polystyrene latex spheres (PSL, specific gravity at 
1.04), suspended in 200-proof ethanol. Based on a first-order estimation, particles ~ 0(1 pm) 
follow the flow well at the applied frequency. The laser was an argon-ion laser operated at 514.5 
nm (green line). The cross-beam half angle between the two transmitted beams was 1.872 
degrees, calibrated at far field (~25 ft) via beam separation measurement. The LDV optics had a 
750.51 mm focal length. The estimated sample volume cross-section was about 175 microns. 
The transmitting probe was set at an angle of 4 degrees down relative to the horizontal plane to 
accommodate getting the sample volume into the slot at 0 mm height without blocking the lower 
beam. The receiving probe was set at an angle of 10 degrees down to ensure that all available 
scattered light was collected unobstructed. 
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The probes were mounted on translation stages to provide lateral and vertical motion. The 
stages have a position accuracy of 1 micron. Once alignment was achieved between the 
transmission and receiving probes, each was moved individually to each new position. Alignment 
was ensured by the signal from the photomultiplier tube. 

LDV signals from the scattering particles were processed using FFT processors. Doppler 
frequency, hence velocity signal, was estimated based on 256 point FFT spectrum analysis at 
each particle arrival. Drive signal waveform was also recorded at the particle arrival time to 
determine the phase information of the velocity signal. 

At each station, 30,000 data points were collected in a period of several minutes, depending 
on the seeding density. Based on the phase information collected at each data point, 30,000 data 
were sorted into 36 phase bins, at 10 degree of phase interval. Within each bin, number of 
samples may vary from a few hundreds to over 2000. Theses samples were used to compute 
phase-locked statistics. 

The uncertainty in the data was calculated using estimates of bias and precision errors in the 
experiment. The estimates were based on a nominal condition using an approximate mean 
velocity in the jet. The bias estimates were based on experimental geometrical parameters, LDV 
processor bias, and biases related to the seeding material used. The total bias and the precision 
were propagated through the coordinate transformation due to the probe downward tilt and 
combined to give an estimate of the total uncertainty in the vertical velocity measured of +/- 0.54 
m/sec. 

Particle Image Velocimetry measurements 

PIV measures an instantaneous velocity field over a grid of points in a plane in the fluid. In this 
test, a Digital PIV system was set up to measure the horizontal and vertical velocity components 
synchronized with the 445 Hz drive signal. The current PIV system includes 1024 x 1280 CCD 
cameras installed with a 200mm Nikon Macro lens. The camera lens was placed approximately 
12” from the laser sheet to cover a field of view about 9 mm wide. The magnification of the 
imaging optics was calibrated with an optical grid target aligned with the laser sheet. The 
accuracy of calibration was within +/-1 pixel over 1240 pixels. Interrogation resolution was set at 
28 by 28 pixels corresponding to about 0.2 mm of measurement area. 

Dual Nd-Yag lasers, operated at 5 Hz rate and 100 mJ output per pulse, were used to 
illuminate a light sheet less than 0.5 mm thick. The laser light sheet was projected perpendicular 
to the slot at the mid-section of the jet slot. Laser pulse separations ( St between double exposure) 
were adjusted between 1.1 ~ 8.0 qsecond at each phase to maximize PIV performance (8 pixels of 
maximum particle displacement). In general, PIV resolves between 1/10 and 1/20 sub-pixel 
resolution. 

Both smoke and polystyrene latex spheres particles were used to seed the entire glass 
enclosure. PIV measurements normally requires higher seeding density than LDV to ensure its 
performance. Smoke particles were produced by a smoke generator using standard smoke fluids 
(specific gravity at 1.022). Smoke particles have a poly-disperse size distribution which may 
cover from sub-micron to tens of micron while PSL particles were are mono-disperse in size 
when manufactured. Although some large agglomerated PSL particles may still be present in the 
flows when sprayed into the test chamber with water-ethanol mixture. Particles larger than 5 pm 
may not follow the jet flows. There was no particle sizer available at the test for in-situ particle 
size distribution measurements. Both particle, smoke and PSL, were tested to determine if any 
effects on the PIV measurements of the jet flows, and there was only slight difference between 
the velocity profiles measured. The glass enclosure was repeatedly seeded to replenish particle 
density for PIV measurements. Four hundred PIV images were taken at 72 phases, in 5 degree 
interval synchronized with the drive signal, to estimate the phase-locked statistics. 
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Flow Field Measurements 


The experiment was conducted separately for the three measurement techniques. A new 
diaphragm was installed in the actuator for each measurement technique due to actuator failures 
or changes in the actuator performance. It should be noted that, the actuator used to generate the 
workshop data (given on the website) was different than the ones used to generate the data given 
in Figs. 3 and 4. The later actuators generated higher jet velocities. However, the overall trends 
are similar. Inside the cavity, the temperature might raise 12-13 degrees F above the ambient air 
after warm-up and then maintained a steady temperature at that point. Hot-wire tests were 
conducted in clean air while the enclosure was moderately seeded for LDV, and heavily seeded 
for PIV test. There was no effect detected of seeding density on cavity pressure waveforms. 

The zero mass jet flow consists of a pulse train of high speed fluids projected form the slot 
accompanied by a pair of vortices. An elongated mean jet emerges when the flow fields are time 
or ensemble averaged. The maximum jet velocity (Fig. 4) was found slightly above the slot, and 
its magnitude may be scaled by the peak cavity pressure, °c (p m a X ) 1/2 - In the near fields, there was a 
close comparison found in the velocity profiles obtained by the LDV and PIV, but their 
magnitudes were off. Good agreements between the LDV and PIV measurements have been 
reported in liquid experiments, but the discrepancy found in this test has not been resolved. On 
the other hand, the hot-wire data showed a faster decay of the jet. In the far fields, however, hot- 
wire and LDV data matched well. Above the maximum jet velocity, the flow is always in upward 
motion, i.e., its height, about 4 slot width, defines the upper limit of the suction action. 



Figure 1: Synthetic jet slot 



Figure 2: Synthetic jet performance with frequency. 


1.1.4 





time 


Figure 3: Waveforms of drive signal (line), cavity pressure (circle), and displacement (square). 
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Figure 4: Maximum, minimum and mean velocity profiles. 
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CASE 1: DETAILS AND SUBMISSION GUIDELINES 


Relevant Details 

Flow issues into an enclosed box (24 inches on each side) with M_freestream = 0. The 
atmospheric conditions varied, but were essentially standard atmospheric conditions at sea level, 
in a temperature-controlled room. These conditions can be given as approximately: p_ambient = 
approx 101325 kg/(m-s A 2) T_ambient inside box = approx 75 deg F (approx 297 K) T inside 
cavity = approx 83 deg F (approx 301 K). Some derived relevant conditions are: density_ambient 
= approx 1.185 kg/m A 3 viscosity_ambient = approx 18.4e-6 kg/(m-s) The diaphragm frequency 
= 444.7 Hz. The 2.0 inch-diameter diaphragm is circular, and is held in place with an o-ring seal 
1.5 inches in diameter. The displacement (D) at the center of the diaphragm as a function of 
phase was provided on the web site. 

The displacement at the center of the diaphragm is offset - it displaces inward LESS and 
outward MORE (i.e., the reference position is not zero) Also in the file, the pressure (P) and 
temperature (T) as a function of phase are given. The location where these quantities are 
measured are as follows. The pressure is measured inside the cavity on the wall opposite the 
center of the diaphragm (the pressure is given with respect to an unspecified reference pressure, 
which may be presumed to be approximately atmospheric). The temperature is measured at the 
bottom (floor) of the cavity. (The voltage input to the diaphragm, prior to being amplified, is also 
included in the file, although this is not relevant information for CFD. It has been offset to 
positive values for LDY purposes.) 

The following image shows a view of the cavity assembly from the underside. The diaphragm 
is on the right side (next to the displacement gauge). The image shows the location of the 
pressure gauge (mounted onto the side wall of the chamber opposite the center of the 
diaphragm), as well as the location of the thermocouple (at the bottom wall of the chamber). 



Synthetic Jet 
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The following figure shows the measured phase-averaged vertical velocity over the center of 
the slot as a function of phase, using three different measurement techniques. The data files can 
be found on the web site. 
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Submission Guidelines 


Numerical predictions of this type of statistically unsteady flow are relatively new. The 
purpose here is to determine the state-of-the-art in modeling these types of unsteady synthetic- 
jet-type flows, so we want to explore which CFD methods work and which do not. The following 
options are available 

■ Either model the internal cavity or specify an unsteady boundary condition at the slot 
exit. 


■ Either model the moving diaphragm, or employ an unsteady boundary condition that 
approximates its effect within the cavity. 
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- Even though the geometry (and flowfield) is inherently three-dimensional, a two- 
dimensional model, across the center of the slot, can be used. 

It is the requirement that contributors detail specifically how they choose to model the case, 
including all boundary conditions and approximations made. As the methodologies used are 
assessed at the workshop, it will be important to know as many details as possible about the 
calculations/simulations. Detailed requirements include: 

1. The case must be run time-accurately, in order to simulate the unsteady nature of the case. 

• la. RANS codes run in time-accurate mode (e.g., URANS) solve directly for phase- 
averaged variables, i.e. <Uj> and cu'iU )> (see Appendix for details). Therefore, RANS 
simulations should result in repeating or very-nearly-repeating periodic solutions. When 
periodicity is achieved, averages over one or more periods of oscillation yields the long- 
time-averaged (time-independent mean) values for these quantities. 

• lb. DNS, LES, or DES simulations will need to be post-processed to obtain both the 
phase-averaged and long-time-averaged (time-independent mean) values. 

2. GRID STUDY : If the case is modeled two-dimensionally, then it is strongly suggested that the 
computation be performed using at least two different grid sizes in order to assess the effect of 
spatial discretization error on the solution. If it is modeled three-dimensionally, then solutions 
using more than one grid size are encouraged, but not required. If more than one grid is used, 
each set of results should be submitted separately. 

3. TIME STEP STUDY: If the case is miodeled two-dimensionally, then it is strongly suggested 
that the computation be performed using at least two different time step sizes in order to assess 
the effect of temporal discretization error on the solution. If it is modeled three-dimensionally, 
then solutions using more than one time step are encouraged, but not required. If more than one 
time step is used, submit each set of results separately. 

Specific quantities that result from the computations at particular locations will be required 
for submission. Note that for all the following, the coordinate system adopted has x across the 
slot (across the 0.05 inch-wide gap) and y up, with the (x,y)=(0,0) origin at the bottom wall of 
the "box" into which the jet issues, at the center of the slot opening. All results from 3-D 
computations are to be given at the z=0 location (across the center of the slot). The requirements 
follow (if it is not possible to provide a particular quantity, simply leave it out of the "variable" 
list, and reduce the number of columns of data submitted): 

a. Long-time-averaged jet width (width is in the x-direction) as a function of vertical height 
above the lower wall, from the wall up to a height of at least 20 mm. Define jet width as the 
horizontal distance between which the vertical velocity exceeds its (maximum+minimum)/2 
value (at that vertical height). Name this file: easel. avgjetwidth.ANYTHING.dat-where 
"ANYTHING" can be any descriptor you choose (should be different for each file if you are 
submitting multiple runs)-the file should be in 2-column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #number of time steps per cycle (pound sign needed) 


1.2.3 



6th line: #brief description of code/method (pound sign needed) 

7th line: #other info about the case, such as spatial accuracy (pound sign needed) 
8th line: #other info about the case, such as turb model (pound sign needed) 

9th line: #other info about the case (pound sign needed) 

10th line: variables="y, mm'V'avg jet width, mm" 

1 1th line: zone t="jet width" 

subsequent lines: y(mm) avgjetw(mm) <- this is the data 
The sample file casel.avgjetwidth.SAMPLE.dat can be downloaded from the web site. 


Example determination of jet width 



x, mm 
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b. Long-time-averaged horizontal velocity (u) and vertical velocity (v) profiles along 

the center of the jet (at the center of the slot), from the wall up to a height of at least 20 mm. 
Also, submit profiles along nine lines parallel to the wall, at heights of 0. 1 mm, 1 mm, 2 mm, 
3mm, 4 mm, 5 mm, 6 mm, 7 mm, and 8 mm above the lower wall, from at least -10 mm to +10 
mm to either side of the jet centerline. Name this file: case l.avgvel.ANYTHING.dat- where 
"ANYTHING" can be any descriptor you choose (should be different for each file if you are 
submitting multiple runs) 

-the file should be in 4-column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #number of time steps per cycle (pound sign needed) 

6th line: #brief description of code/method (pound sign needed) 

7th line: #other info about the case, such as spatial accuracy (pound sign needed) 

8th line: #other info about the case, such as turb model (pound sign needed) 

9th line: #other info about the case (pound sign needed) 

10th line: variables="x, mm","y, mm","u, m/s","v, m/s" 

11th line: zone t="centerline" 

subsequent lines: x(irnn) y(mm) u(m/s) v(m/s) <- this is the data along x=centerline 
next line: zone t="y=0. 1 mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y=0. 1mm 
next line: zone t="y=l mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y= 1mm 
next line: zone t="y=2 mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y=2mm 
next line: zone t="y=3 mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y=3mm 
next line: zone t="y=4 mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y=4mm 
next line: zone t="y=5 mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y=5mm 
next line: zone t="y=6 mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y=6mm 
next line: zone t="y=7 mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y=7mm 
next line: zone t="y=8 mm" 

subsequent lines: x(mm) y(mm) u(m/s) v(m/s) <- this is the data along y=8mm 
The sample data file casel.avgvel.SAMPLE.dat can be downloaded from the website. 

c. Phase- averaged quantities at 8 different phases during the cycle: 0 deg, 45 deg, 90 deg, 135 
deg, 180 deg, 225 deg, 270 deg, 315 deg.; where you should align the phases of your 
computation as described below. Submit the following phase-averaged <> quantities: u (m/s), v 
(m/s), u'u'bar (m A 2/s A 2), v'v'bar (m A 2/s A 2), and u'v'bar (m A 2/s A 2), where u = phase-averaged 
horizontal velocity component (parallel to lower wall, across the slot), v = phase-averaged 
vertical velocity component (up, away from lower wall) u'u'bar = phase-averaged turbulent 
normal stress in horizontal direction (optional) v'v'bar = phase- averaged turbulent normal stress 
in vertical direction (optional) u'v'bar = phase- averaged turbulent shear stress in x-y plan The 
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locations for these data are the same as for the long-time-averagedquantities. Name these files: 
casel.phaseOOO.ANYTHING.dat 

case 1 ,phase045 .ANYTHING.dat 
case 1 .phase090.ANYTHING.dat 
case 1 .phase 1 35 .ANYTHING.dat 
case 1 .phase 1 80.ANYTHING.dat 
case 1 ,phase225 .ANYTHING.dat 
easel ,phase270. ANYTHING.dat 
case 1 ,phase3 1 5 .ANYTHING.dat 

-where "ANYTHING" can be any descriptor you choose (should be different for each file if you 
are submitting multipleruns) 

-the file should be in 7-column format: 


1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #number of time steps per cycle (pound sign needed) 

6th line: #brief description of code/method (pound sign needed) 

7th line: #other info about the case, such as spatial accuracy (pound sign needed) 
8th line: #other info about the case, such as turb model (pound sign needed) 

9th line: #other info about the case (pound sign needed) 


variables: 


m/s","v, m/s". 


m A 2/s A 2","vv, 


m A 2/s A 2","uv,m A 2/s A 2" 


11th line: zone t="centerline" 


subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along x=centerline 
next line: zone t="y=(). 1 mm" 
subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=0. 1mm 
next line: zone t="y=1 mm" 
subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=lmm 
next line: zone t="y=2 mm" 
subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=2mm 
next line: zone t="y=3 mm" 
subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=3mm 
next line: zone t="y=4 mm" 
subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=4mm 
next line: zone t="y=5 mm" 
subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=5mm 
next line: zone t="y=6 mm" 
subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=6mm 
next line: zone t="y=7 mm" 
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subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=7mm 
next line: zone t="y=8 mm" 
subsequent lines: x(mm) y(mm) u(m/s) v(m/s) 
uu(m A 2/s A 2) vv(m A 2/s A 2) uv(m A 2/s A 2) <- this is the data along y=8mm 

The sample datafile casel.phaseOOO.SAMPLE.dat can be downloaded from the website. 


d. Phase-averaged time-history values of <u> and <v> as a function of phase (deg) at three 
approximate point locations: (x,y)=(0mm,0.1mm), (Omm, 2mm), and (1mm, 2mm). Give the data 
at every time step taken. In other words, if your time step yields 100 steps per cycle, then give 
100 phases between 0 deg and 360 deg. You should align the phases of your computation as 
described below. Name this file: casel.phasehist.ANYTHING.dat -where "ANYTHING" can be 
any descriptor you choose (should be different for each file if submitting multiple runs) 

-the file should be in 5 -column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #number of time steps per cycle (pound sign needed) 

6th line: #brief description of code/method (pound sign needed) 

7th line: #other info about the case, such as spatial accuracy (pound sign needed) 

8th line: #other info about the case, such as turb model (pound sign needed) 

9th line: #other info about the case (pound sign needed) 

10th line: variables="phase, deg","x, mm","y, mm","u, m/s","v, m/s" 

11th line: zone t="x=0 mm, y=(). 1 mm" 

subsequent lines: phase x(mm) y(mm) u(m/s) v(m/s) <- this is the data at x=0mm, 
y=0.1mm 

next line: zone t="x=0 mm, y=2 mm" 

subsequent lines: phase x(mm) y(mm) u(m/s) v(m/s) <- this is the data at x=0mm, y=2mm 
next line: zone t="x=l mm, y=2 mm" 

subsequent lines: phase x(mm) y(mm) u(m/s) v(m/s) <- this is the data at x=lmm, y=2mm 
The sample datafile casel.phasehist.SAMPLE.dat can be downloaded from the website. 


e. Field line-contour-plots (in one of the following formats: ps, eps, or jpg) of phase-averaged 
<u>-velocity and <v>-velocity at the following phases: 45 deg, 90 deg, 135 deg, and 225 deg.; 
where you should align the phases of your computation as described below. These plots should 
be black-and-white line plots only. The plots should go from approx x=-3mm to 3 mm, and 
y=0mm to y=8mm. The x-to-y ratio of the plot should be 1.0. For all plot files, plot line 
contours of -20 m/s through 30 m/s in increments of 1 m/s. Either (a) plot the lines of negative 
velocity as dashed lines and the lines of positive velocity as solid lines, or (b) label the contour 
lines, or (c) do both. The purpose of submitting these plots is to get a qualitative picture of the 
phase-averaged flowfield at particular selected times of interest. Altogether, submit 8 plots files. 
Name these files: easel. phase045.u. ANYTHING.eps 
easel .phase090.u. ANYTHING.eps 
casel.phasel35.u-ANYTHING.eps 
case 1 .phase225 .u. ANYTHING.eps 
case l.phase045. v.ANYTHING.eps 
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easel .phase090. v.ANYTHING.eps 
case 1 .phase 1 3 5 . v . AN YTHING. eps 
casel.phase225.v-ANYTHING.eps 

(where the "eps" in this case means encapsulated postscript - use ps, or jpg instead if 
appropriate). 

Definition of Phase for the Computations 

Matching the same phases with the experiment is not necessarily straightforward. One way to do 
it is to try to align a quantity from experiment (such as vertical velocity near the jet exit), but this 
can be imprecise because the CFD and experimental data are not necessarily well-behaved sine- 
waves. The best criteria for determining phase may in fact be a different measure altogether. 

On the other hand, in this workshop a goal is to be able to compare CFD results with each other, 
so it is important to try to achieve the same phase definitions in order that all computations are 
similarly aligned. We have tried to choose criteria for determining phase that approximates 
experiment AND is specific enough so that different CFD solutions can be meaningfully 
compared. 

Therefore, although participants are given some latitude to determine phase as appropriate, we 
encourage everyone to use the following steps to define phase in a uniform fashion for Case 1 : 

• Step 1. Output vertical phase-averaged velocity (v) at the following point in space (over 
the slot) as a function of your iteration or time step number: (x,y,z)=(0, 0.1, 0)mm. Find 
vmax and vmin over the course of one phase-averaged period. 

• Step 2. Compute the mid-value vavg=(vmax+vmin)/2 

• Step 3. Define Phase=340 deg as the time when your velocity at this location 
approximately equals vavg (INCREASING). All other phases can be referenced from 
this, via the following relation: 

Phase=(iter-it340)*360/nstep+340 

where: 

iter iteration (or time step) number 
nstep=no. of time steps per cycle 

it340=iteration number when Phase=340 according to the above criteria 

For example, if the calculation is running at 360 steps per cycle, in order to match the above 
criteria at time step number 5575, then 

Phase=iter-5235 

Thus, when iter=5235, then Phase=0; when iter=5415, then Phase=180; when iter=5595, then 
Phase=360. Note that Phase=360 also corresponds with Phase=0 (it repeats every 360 deg). This 
is illustrated in the following figure: 
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Example of defining phase from phase -averaged 



As another example, if the calculation is running at 1080 steps per cycle, in order to match the 
above criteria at time step number 10002, then 

Phase=(iter- 1 0002)/3+340 

Thus, when iter=8982, then Phase=0; when iter=9522, then Phase=180; when iter=10062, then 
Phase=360. 
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APPENDIX: FLOWS WITH IMPOSED PERIODIC FORCING 


Mathematical Background 


The flow field is induced by an oscillatory jet of air emanating from a solid boundary. This oscil- 
latory motion introduces a wave pattern that can be used as a timer for selective sampling of the 
flow field fluctuations. Thus, the weak organized wave motion can be extracted from a background 
field of finite turbulent fluctuations. The mathematical framework outlined here follows closely that 
originally presented by Hussain and Reynolds (1970, 1972). Some previous work in this area that 
may also be of interest are reported in Gatski and Liu (1980) and Obi, Ishibashi and Masuda (1997). 

If such traveling waves exist within the fluctuating flow field, any fluctuating quantity /(x, t) 
can be decomposed into 

/ Cm) = / (x) + / (x,i) + /' (x,t) , (i) 

where / (x) is the (time) mean value, / (x, t) is the statistical contribution of the organized motion, 
and // (x, t) the turbulence. Since the flow is not stationary, the ergodic hypothesis does not hold 
and the ensemble mean does not equal the long time average. Thus, a time average is defined as 


/» 


lim — 

X — »oo T 


f (x, t)dt 


( 2 ) 


and a phase average is defined as 


1 N 

(/ ( x ^)) = J J™ 0 ]v^ /(x ’ i + rar) 

n= 0 

(3) 

where r is the period of the imposed oscillatory motion. The phase average is the average at any 
point in space of the values of / that are realized at a particular phase cf) in the cycle of the oscillating 
jet. The wave component is then given by 

/ (x,£) = (/ (x,i)) - / (x) , 

(4) 

so that 

/(x, t) = (/(x,t)) + /' (x,t), 

(5) 


These results show that at any point in the flow a time varying function /(x, t) can be partitioned 
into the three components parts defined in Eqs. (2), (3), and (4) given a reference oscillating signal 
at a given frequency (period). The following properties are associated with the triple-decomposition 
and prove useful in extracting information about both the mean field(s) and the statistical correla- 
tions of the fluctuation field. These relations include 

fg = fg , (fg) = f(g), < fg > = / (g) , if) = (/} = f, 

— — — / ~ r (6) 

/ = 0, /' = 0, </'> = 0, /y=(/y) = 0. 

These relations show that the fluctuation f f is centered about both the time average and phase 
averaged means and that the background turbulence is uncorrelated with the organized motion. 
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Phase-Averaged Navier-Stokes Equations 

In general, both pressure-velocity and density-velocity numerical solvers are used in RANS CFD 
codes. In this section, the phase-averaged formulations for both the incompressible and compress- 
ible forms of the Navier-Stokes equations will be presented. For contributors using time-accurate 
RANS-type formulations, the dependent variables calculated are the phase-averaged variables de- 
fined above and shown below to satisfy transport equations that are formaly equivalent to the usual 
RANS equations. 

Incompressible Phase-Averaged Equations 

With the decomposition and mathematical framework outlined in Section , it is possible to develop 
a set of phase-averaged Navier-Stokes equations that describe the behavior of the phase-averaged 
mean quantities as well as the corresponding equations for the phase-averaged turbulent quantities 
(Gatski and Liu, 1980). 

In Section , the phase-averaged mean equations are presented and in Section the phase-averaged 
turbulent quantities are presented. 

Phase-averaged mean equations 
The mean flow equations are given by 

9 ( m ») = n 

dx t 

D (3) = 9 ( u i) | / \ d ( u i) = 1 d(p) d 2 (uj) d ( u i u 'j) 

Dt dt 3 dx'0' p dxi dx 2 dx 3 

From these results, it is possible to obtain the time-independent mean values, for example 

Ui = ( Ui ), w'w' = (u' t u'P (8) 


(7a) 

(7b) 


Phase-averaged turbulent equations 


As an example, the phase-averaged formulation for a two-equation (K)-(e) model will be given. 
Other models under consideration can be formulated in terms of phase-averaged variables in a sim- 
ilar fashion. 

The transport equations for the phase-averaged turbulent kinetic energy (K) (= (n'n') /2) and 
dissipation rate (e) are given by 


D(K) 

Dt 


[ u Wk\ 


d_(ui 

dxb 


-(e) + 


d_ 

dxi 


, Vt 

v H 

ctk 


dJKY 

dxb 


( 9 ) 


D(e) 


= -c £l 


(K) 


[U:U h 


d{ui 

dxi. 


-c, 


£ 2 


+ 


d 


(K) dx k 


, Vt 

v H 

<J £ 


d(s) 


dxi 


( 10 ) 


Dt 

where ax, cr £ , ( ' \ . and ( ' ■> are closure coefficients, and an eddy viscosity relationship is assumed 
for the individual phase-averaged stress components so that 


with 




and CV a closure coefficient. 


3 (K)^j M q x \ ! 


+ 


3x; 


v t = C l 


(i<y 


(Ha) 

(lib) 
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Compressible Phase-Averaged Equations 

In a compressible formulation, where Favre averaged variables are used, that is 'pul/ p, the definition 
of the phase average has to be altered. In this case the compressible phase-averaged variables are 
defined as 

\ u i) — > (12) 

and the flow variables are partitioned as in Eq. (5). In Section , the phase-averaged mean equations 
are presented and in Section the phase-averaged turbulent quantities are presented. 

Phase-averaged mean equations 

In terms of phase-averaged Favre variables, the mass conservation equation can be written as 


D(p) _ d(p) 


d ip) = _ / \ 9 (uj) 

dx 3 w dx 3 


and the conservation of momentum as 


D(ui) _ dip) d {oij) 9 (p)( u i u j ) 


_| 

dxj dx;. 


where 


/ \ OT7 f 1 ( d i u i) , d i u 3 ) 

(< ’« >=2 '‘ 2(-&r + “te 


1 d(u k ) ^ 
3 dxb 


A phase-averaged conservation of total energy equation can be written as 

9({p) m d{ ( Uj ) (p) (H)) _ 0 (upE) d (qj) 
dt dx 3 dx 3 dx 3 5 

where the total energy and total enthalpy, (p) ( E ) and (p) (H), respectively, are 

(£) = 

(H) = (E) + M, 


in) = - { k T -K— ) - - k T 


r d i T ) 


(ujpE) = c p (p) (u'jT') + (u t ) {{p) (u'iu'j) - {(Tij)) + 


ip) (<'<'< 


\ U i <J ij / 


From these results, it is possible to obtain the time-independent mean values, for example 


p = (/?), u t = (m), E = (E), q % = (q t ) % t/'n' = 
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Phase-averaged turbulent equations 


Once again, a phase-averaged formulation for a two-equation (K)-(e) model is used. Other models 
can be formulated in terms of phase-averaged Favre variables in a similar fashion. 

The transport equations for the phase-averaged turbulent kinetic energy (K) (= (n'n') /2) and 
dissipation rate (e) are given by 


(p) 


D(K) 

Dt 


( P ) < u 'i u 'k ) ~ (p) ( e ) + 


dxi 


P> + 


Pt \ d(KY 


(JK 


dxi 


(18) 


/ \ D (e) (e) 

w ~dt = ~ c -' « 


z , , \ d ( Ui ) . . (e) 

\ U i U k) ~^Z C e2 ( P ) 7777 + 


d_ 

dxi 


fl + 


fi t \ d(s) 


d xi 


(19) 


(K) ' 1 k/ dx k ez ^'(K) 

where a k, <j £ , C £ i , and ( V 2 are closure coefficients, and an eddy viscosity relationship is assumed 
for the individual phase-averaged stress components so that 


(P) <“'<'> = \ (P) (K) Sij - fi t 


+ 


dx; 


(20a) 


with 

(K) 2 

Pt = CV (p) Vt- (20b) 

\ € ) 

and C M a closure coefficient. 
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CASE 1: NUMERICAL SIMULATION OF A SYNTHETIC JET INTO 

A QUIESCENT AIR 


I. Mary 1 


1 CFD and Aeroacoustic Departement, ONERA Chatillon, France 


Introduction 

Three different simulations of the case 1 have been carried out in order to assess the influence of the tur- 
bulence modelling. Firstly two-dimensional URANS and laminar simulations have been first investigated. 
As the Reynolds number based on cavity charateristics (like the diaphragm velocity and the exit size of the 
cavity) is quite small (less than 10000), it is not obvious that an URANS simulation will be better adapted 
to this flow case than a laminar simulation. Moreover the discontinuities in the geometry of the cavity could 
create complex flow patterns inside the cavity despite the low value of the Reynolds number. Thus three- 
dimensional Large Eddy Simulation have been realized for the same two-dimensional geometry in order to 
evaluate the intensity of the three-dimensional effects inside the cavity. 


Solution Methodology 

General description 

The solver FLU3M, developed by ONERA, is based on a cell-centered finite volume technique and struc- 
tured multi-block meshes. The viscous fluxes are discretized by a second-order accurate centered scheme. 
For efficiency reason, an implicit time integration is employed to deal with the very small grid size en- 
countered near the wall. Then a three-level backward differentiation formula is used to approximate the 
temporal derivative associated with the unsteady Navier-Stokes equations, leading to second-order accu- 
racy. An approximate Newton method is employed to solve the non-linear problem. At each iteration of this 
inner process, the inversion of the linear system relies on Lower-Upper Symmetric Gauss-Seidel (LU-SGS) 
implicit method. More details about these numerical points are available in Ref.[l]. 

Euler fluxes discretization 

Usually LES requires a high-order centered scheme for the Euler fluxes discretization (with spectral-like 
resolution) in order to minimize dispersive and dissipative numerical errors. However such scheme cannot 
be applied easily in complex geometry. Indeed, most of aerodynamic codes able to deal with such a geometry 
are based on finite-volume technique in order to handle degenerated cell. Thus getting a high-order method 
becomes very time consuming due to the high-order quadrature needed to compute the fluxes along the cell 
boundaries. As several works (see Ref. [2] for instance) have shown that LES can be carried out with low- 
order centered scheme in case of sufficient mesh resolution, only second-order accurate scheme is employed 
in this study. But a special effort has been carried out to minimize the intrinsic dissipation of the scheme 
and its computational cost. 

The AUSM+(P) [3, 4] scheme, whose dissipation is proportional to the local fluid velocity, constitutes 
the basis of the discretization, because it is well adapted to low Mach number boundary layer simulations. 
However several modifications have been introduced to enhance its accuracy and computational cost. As we 
are not interested by the shock capturing properties of the scheme, simplified formula are developed for this 
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study to approximate the Euler fluxes: 

Fl +1/ = Ui(Q l + Q r )/ 2 - \U dis \(Q R - Q l )/ 2 + P ( 1 ) 

where U\ denotes the interface fluid velocity, L/R the left and right third-order MUSCL interpolation. 
The state vector Q is defined as (p, pu\, puz, pu?,, pE + p )*, whereas the pressure term P is given by 
(0, (Pl+Pr)/ 2, 0, 0, 0). The symbol (7^ s , which indicates a local fluid velocity, characterizes the numerical 
dissipation acting on the velocity components. To enforce the pressure/velocity coupling in low Mach 
number zone, a pressure stabilization term is added to the interface fluid velocity as done by Rhie and Chow 
[6] for incompressible flow: 

Ui = (u 1L + uir)/2 - c 2 (pR -Pl) (2) 

where C 2 is a constant parameter. The parameter Udi s is defined as follow: 

U dis = $ max {\ui L + Ui R \/2, ci) (3) 


where c\ is a constant parameter and a sensor used to minimize the numerical dissipation. For an accuracy 
reason, the values of c\ and C 2 should be chosen as small as possible to minimize the numerical dissipation. 
For a stability reason, these parameters cannot be smaller than the threshold value 0.04, which has been 
determined in [5]. The present sensor $ is a binary function, which only depends on the smoothness of the 
primitive variables ^ = (p, fq, U 2 , u^^p) 1 . If no spurious oscillation is detected on i/; at cell i, is equal 
to 0. Otherwise 4> is set to 1 at % + 1/2 and i — 1/2 interfaces. To detect a spurious oscillation and define 
exactly the sensor used in Eq.(3), the following functions are introduced: 


A i f-i if (< 4 + 2 - <; k+i)(<i>i+i ~4>i) < o 

v \ 1 else 

W» 4 1 if A k+ <°” A k+ A £‘< o 

Wk 1 0 else 

$ = max{Wp k ) for k= 1 ,5 


(4) 


For the two-dimensional laminar and URANS simulations, there is no need to capture the small scales 
eddies. Therefore for these simulations the value of $ is fixed to 1, and the scheme is then purely upwind. 


Model Description 

Three different kinds of turbulence modelling have been used to compute the test case 1 in order to as- 
sess their influence for such flow. First a laminar two-dimensional computation has been carried out. A 
two-dimensional URANS simulation has been performed as well, based on the classical Spalart-Allmaras 
model. Finally a three-dimensional FES simulation has been computed. This simulation is based on an eddy 
viscosity model developed at ONERA: the selective mixed scale model. 


Selective mixed scale model 

Due to the large computational cost of the intended simulation, only one explicit subgrid scale model has 
been used for this study. The selective mixed scale model, developed by Sagaut and Fenormand, has been 
retained, because a good compromise between accuracy, stability and computational cost is obtained [7]. 
More particularly, the use of a selective function allows to handle transitional flows [8], which is one of the 
keypoint of the present application. The eddy viscosity is given using a non-linear combination of the filtered 
shear stress tensor, a characteristic length scale, the small scale kinetic energy, and a selective function: 

/it = pfdo (0)Cm A 3/2 ^0-5 S~j (u)Sij (u) VK: (5) 
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Figure 1 : Diaphragm and grid representation 


where the test field kinetic energy is evaluated as Kc — y0.5 (ui)' (ui)' . The test field (ui)' is extracted 
from the resolved field, (u t )' = ui — Ui, by employing a three-dimensional averaging test filter denoted by 
Ui = Ai(Aj(Ak(ui))), with the following definition of A\\ 

71/(0) = O.250/_i + 0.50/ + O.250/+1 (6) 


The characteristic length scale, A, is given by the cube root of the cell volume of the mesh, and the constant 
parameter, C m = 0.06, is derived in Ref. [9] from isotropic homogeneous turbulence flow cases. The 
selective function is defined by: 

f (f))= n ii 9 > 9o = 10° 

| fan 4 (|)/fan 4 (^-) otherwise 

where 6 is the angle between the filtered vorticity, a; (it) and the local averaged filtered vorticity, 


Implementation and Case Specific Details 

Diaphragm modelling 

A two-dimensional geometry of the case 1 has been retained for this study (see the Fig. 1). No motion is 
applied to the mesh at the interface of the diaphragm. Therefore the velocity is imposed at the boundary 
of the computational domain by using a time derivative of the time history of the diaphragm displacement 
given in the workshop experimental data. 

V(y) — 0.5sm(2 x n x 444.7 x t) My (8) 

The use of this unsteady boundary condition generates vortices inside the cavity, which are illustrated on the 
Fig. 2. 

Concerning the thermodynamic variables, the density and the pressure are obtained at the wall thanks to 
a zero order extrapolation. 
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Figure 2: Illustration of the flow structures near the diaphragm: a) t=T/4. b) t=T/2. c) t=T/3. d) t=T. 


Time integration 

The time step of all the simulations is fixed at At = 0.4497414/is (5000 time steps per period of the 
diaphragm motion). The number of iteration in the Newton process is equal to 4, which is sufficient for this 
case to deacrease the residual of the non-linear equations by about one order of magnitude. 

Grids 

A multi-block mesh has been generated for this study. The size of the two-dimensional computational 
domain is identical to those provided on the web-site of the workshop. However the cell distribution differs 
significantly. At the cavity exit, the grid spacing is equal to 0.1 mm in the y direction. In the exit channel 
of the cavity, there are 80 cells in the x direction. The grid spacing at the wall is equal to Ax = 0.004 
mm. Therefore there are around 30 cells located in the boundary layer. In the streamwise direction, there 
is 134 cells distributed in the upper part of the channel (see Fig. 1). For the external flow there are 217 and 
131 cells distributed in the x and y direction, respectively. A refinement has been used near the cavity exit. 
Therefore in the y direction, 80 cells are located between y=0 and y=20 mm, whereas in the x direction the 
mesh has been adapted to the tickening of the jet. Therefore 100 cells are located inside the jet diameter 
independantly of the y location. For the three-dimensional LES simulation, the same grid has been used in 
the (x,y) plan. In the z direction, the span wise extent of the computational domain is equal to 1^ = 1.3mm, 
which is almost the size of the jet diameter at the cavit exit. In the spanwise direction 20 cells are uniformely 
distributed. 

Non-reflective boundary conditions are used in the farfield. At the wall, a no-slip adiabatic condition 


1.3.4 



has been retained. For the three-dimensionnal LES simulation, the same boundary conditions are used. 
Moreover flow periodicity is assumed in the spanwise direction. 
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Governing Equations and Solution Methodology 

In the present work, a generalized form of the thin-layer Navier-Stokes equations is used to model the 
turbulent, viscous flow encountered in a synthetic jet. The equation set is obtained from the complete Navier- 
Stokes equations by retaining the viscous diffusion terms normal to the solid surfaces in every coordinate 
direction. For a body-fitted coordinate system (£, 77 , Q fixed in time, these equations can be written in the 
conservative form as: 

T fl(U) , d(F - F v ) (9(G - G v ) fl(H-H v ) _ 

dt + + dr] + d( U 

where U represents the conserved variable vector. The vectors F, G, H, and F v , G v , H v represent the 
convective and diffusive fluxes in the three transformed coordinate directions, respectively. In Eqn. (1), J 
represents the cell-volume or the Jacobian of coordinate transformation. A multigrid based, general purpose 
multi-block structured grid approach is used for the solution of the governing equations. In particular, 
the TLNS3D flow code is used in this study to solve Eqn. (1). Several references exist describing the 
discretization and algorithmic details of TLNS3D code. We include only a brief summary of the general 
features, and refer to the work of Vatsa and co-workers [1,2] for further details regarding the TLNS3D code. 


Numerical Discretization 

The spatial terms in Eqn. (1) are discretized using a standard cell-centered finite- volume scheme. The con- 
vection terms are discretized using second-order central differences with scalar/matrix artificial dissipation 
(second- and fourth- difference dissipation) added to suppress the odd-even decoupling and oscillations in 
the vicinity of shock waves and stagnation points [3,4]. The viscous terms are discretized with second-order 
accurate central difference formulas [1]. For the present computations, the Spalart-Allmaras (SA) model [5] 
and the Menter’s SST model [ 6 ] are used for simulating turbulent flows. 

For temporal discretization, the convective and diffusive terms are grouped together, and Eqn. (1) can 
be rewritten as: 


— = -C(U) + D p ( U) + D a ( U) (2) 

where C(U), D p ( U), and D a ( U) are the convection, physical diffusion, and artificial diffusion terms, 
respectively. The cell-volume or the Jacobian of coordinate transformation being included in these terms. 
The time-derivative term can be approximated to any desired order of accuracy by the Taylor series 

— a iU n + (I 2 U n 1 + ci^U n 2 + ...] (3) 
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Eqn. (3) represents a generalized backward difference scheme (BDF) in time, where the order of accuracy is 
determined by the choice of coefficients a 0 , Ui, a 2 ... etc. For example, if Ho = l.|| a\ — —2 and a 2 = .5, 
it results in a second-order accurate scheme (BDF2) in time, which is the primary scheme chosen for this 
work due to its robustness and stability properties [7]. 

In Eqn. (3), the superscript n represents the time level at which the solution is known, and n+1 refers 
to the next time level to which the solution will be advanced. Similarly, n- 1 refers to the solution at one 
time level before the current solution. Regrouping the time-dependent terms and the original steady-state 
operator leads to the equation: 


^ un+ i , 

At At 


S( U n ) 


( 4 ) 


where £’(U n,n_1 ’-) depends only at solution vector at time levels n and prior, and S represents the steady 
state operator or the right hand side of Eqn.(2). By adding a pseudo-time term, we can rewrite the above 
equation as: 


ao +1 

dr At At 


S( U n ) 


( 5 ) 


Solution Algorithm 

The algorithm used here for solving unsteady flows relies heavily on the steady- state algorithm available 
in the TLNS3D code [1]. The basic algorithm consists of a five-stage Runge-Kutta time-stepping scheme 
for advancing the solution in pseudo-time, until the solution converges to a steady-state. Efficiency of this 
algorithm is enhanced through the use of local time-stepping, residual smoothing and multigrid techniques 
developed for solving steady-state equations. In order to solve the time-dependent Navier-Stokes equations 
(Eqn. 5), we added another iteration loop in physical time outside the pseudo-time iteration loop in TLNS3D. 
For fixed values of f?(U n,n_1 ’"), we iterate on U n+1 using the standard multigrid procedure of 

TLNS3D developed for steady-state, until desired level of convergence is achieved. This strategy, originally 
proposed by Jameson [8] for Euler equations and adapted for the TLNS3D viscous flow solver by Melson 
et. al [7] is popularly known as the dual time-stepping scheme for solving unsteady flows. The process is 
repeated until the desired number of time-steps have been completed. The details of the TLNS3D flow code 
for solving unsteady flows are available in Refs. [7] and [9]. 


Boundary Conditions 

Except for the moving diaphragm, standard boundary conditions, such as the no-slip, no injection, fixed 
wall temperature or adiabatic wall, far-held and extrapolation conditions are used as appropriate for the 
various boundaries. The most accurate procedure to simulate the moving diaphragm would require moving 
grid capability. For simplicity, we chose to simulate this type of boundary condition by a periodic tran- 
spiration velocity condition. The frequency of the transpiration velocity at the diaphragm surface in the 
numerical simulation corresponded to the frequency of the oscillating diaphragm, and the peak velocity at 
the diaphragm surface was obtained from numerical iteration so as to match the experimentally measured 
peak velocity of the synthetic jet emanating from the slot exit. The pressure at the moving diaphragm is also 
required for closure. However, in the absence of unsteady pressure data from the experiment, we imposed 
a zero pressure gradient at the diaphragm boundary. We also tested the pressure gradient boundary condi- 
tion obtained from one-dimensional normal momentum equation [10], which had very little impact on the 
solutions. Due to the simplicity and robustness, we selected the zero pressure gradient boundary condition 
at the diaphragm surface. 
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Implementation and Case Specific Details 

Although the actuator geometry is highly three-dimensional, the outer flow field is nominally two-dimensional 
because of the high aspect ratio of the rectangular slot. For the present study, this configuration is modeled 
as a two-dimensional problem. A multi-block structured grid available at the CFDVAL2004 website is used 
as a baseline grid. 

The periodic motion of the diaphragm is simulated by specifying a sinusoidal velocity at the diaphragm 
surface with a frequency of 450 Hz., corresponding to the experimental setup. The amplitude is chosen so 
that the maximum Mach number at the jet exit is approximately 0.1, to replicate the experimental conditions. 
At the solid walls zero slip, zero injection, adiabatic temperature and zero pressure gradient conditions are 
imposed. In the external region, symmetry conditions are imposed on the side (vertical) boundaries and 
far-field conditions are imposed on the top boundary. A nominal free-stream Mach number of 0.001 is 
imposed in the free stream to simulate incompressible flow conditions in the TLNS3D code, which solves 
compressible flow equations. 

The code was run in unsteady (URANS) mode until the periodicity was established. The time-mean 
quantities were obtained by running the code for at least another 15 periods and averaging the flow quantities 
over these periods. The phase-locked average of flow quantities were assumed to be coincident with their 
values during the last full time period. 


Results/Discussion 

Only a brief set of results is included here, detailed results will be available in the workshop proceedings. 
Majority of the solutions were obtained with SA turbulence model on the baseline grid with a time- step of 
At = T per i 0 ( t/ 72, where T per i 0 d represents the physical time for a complete period. In addition, results 
were generated for the following conditions: (a) coarser grid (eg) obtained by eliminating every other point 
from the baseline grid; (b) finer grid (fg) generated by increasing resolution in the normal direction by 50 
(c) lower time-step (lower dt), where At = T per i 0 d/ 144 ; and (d) Menter’s SST turbulence model. The 
results for the vertical velocity near the slot exit center line (x=0., y=0.1 mm) from these computations are 
compared in Fig. 1, which shows very little sensitivity to turbulence model change, and spatial and temporal 
resolutions. 

Next we show the results for time-averaged velocity along the jet center line in Fig. 2. In the proximity 
of the slot exit (y < 2 mm ), all 5 sets of results are indistinguishable, and are in good agreement with the 
hotwire data. However, further away from the slot, the coarse grid (eg) results show significant deviation 
from the baseline results. The effect of turbulence model type (SA or SST) are also more pronounced in the 
outer region. Similar trends are observed in Fig. 3 for the jet width based on time-averaged solutions. 

Based on these results and more detailed comparisons (not shown here) of the computational results 
for time-averaged and phase-locked average velocities, the effect of lowering the time-step and refining the 
normal mesh were found to be very small. Noticeable differences were observed for the coarse grid (eg) 
results in the outer regions. The difference between the SA and SST turbulence model results become more 
pronounced in regions further away from the slot exit. 
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Figure 1: Time- variation of V- velocity at slot exit, x=0., y=0.1 mm 
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Figure 2: Center-line average velocity comparisons 



Figure 3: Jet-width comparisons 
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Introduction 

Computations of a synthetic jet (Case 1) are performed with Code.Saturne , a finite-volume solver on un- 
structured grids developed at EDF [1]. The purpose of the present work is to investigate the ability of stan- 
dard turbulence models (eddy-viscosity and Reynolds-stress models) to close the phase-averaged Navier- 
Stokes equations. Therefore, this unsteady flow was run with usual RANS-equations solved in time-accurate 
mode (URANS), using the standard k-e model [4] and the Rotta+IP second moment closure [3, 2]. Wall 
functions were used at all solid boundaries. All the computations are 2-dimensional. 

In order to investigate the influence of the grid and the time step, three runs were performed with the 
standard k-e model and one with the second moment closure. A summary of the characteristics of these 
computations is given in the table below. Note that the time steps are given in degrees, i.e., the value 
corresponds to 360/ At, where / is the frequency of the jet. 



Model 

Convection scheme 

Time marching 

grid 

time step (deg) 

1 

Standard k-e 

2nd order 

2nd order 

coarse grid 

0.5 

2 

Standard k-e 

2nd order 

2nd order 

coarse grid 

0.25 

3 

Standard k-e 

2nd order 

2nd order 

fine grid 

0.25 

4 

RSM 

2nd order 

1st order 

coarse grid 

0.5 

5 

RSM 

2nd order 

1st order 

coarse grid 

0.125 


Solution Methodology 

The 2-D computational domain is of the same extent as the box in which the oscillating jet emanates in the 
experiment in order to avoid undesired confinement effects. The internal cavity that generates the flow is 
not modeled: experimental values are imposed as unsteady boundary conditions at the slot exit. 

Computations are started from rest state with small values for turbulence quantities and run until very- 
nearly-repeating periodic solutions are obtained: periodicity is satisfactorily achieved after about 10 cycles. 
Note that periodicity was only evaluated in the region where the flow is of interest (basically the region 
in which the results can be compared with experiments): indeed, the domain is so large that it would take 
hundreds of cycles to obtain a periodic solution in the entire domain. 

Computations with the k-e model were computed with second-order schemes (central differencing, 
Adams-Bashford): even if the computations were started abruptly form rest state, with a coarse time step 
for the first cycles (360/ At ~ 5, i.e., CFL ~ 10), it was not necessary to switch to first-order schemes, 
which shows the robustness of the code. However, with the second moment closure, the computations had 
to be started using upwind differencing and first order time marching, and restarted after 5 cycles with central 
differencing. Due to a shortage of time, results using second-order time marching have not been obtained: 
however, the time step was reduced in computation 5 to limit numerical error. 
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Figure 1: Coarse grid. Left: view of the whole domain. Right: view of the slot exit region. In both views, 
the slot is materialized by two vertical lines. 


Implementation and Case Specific Details 

Boundary conditions 

• Standard wall functions are used on the solid walls of the box. 

• Unsteady Dirichlet boundary conditions are imposed at the slot exit: 

- phase-averaged vertical velocity over the center of the slot as a function of phase is taken from 
the LDV experimental data (at y = 0.1mm). Note that the velocity is uniform across the slot 
exit. 

- phase-averaged v' over the center of the slot as a function of phase is taken from the hot wire 
experimental data (x = 0, y — 0 , 2 = 0). For the k-e model, k = 3/2 i/ 2 is imposed. For the 
Reynolds-stress model, u 2 = v 2 = w 2 = v’ 2 and uv = 0 are imposed. The inlet dissipation rate 
is calculated at every time step such that the turbulent integral length scale equals half the slot 
width. 

• Since the flow is considered incompressible, an outlet is necessary to allow mass conservation at each 
time step. Two small outlets are specified at the lower corners of the box. When the jet is blowing, 
these regions correspond to standard outlet conditions. During the aspiration phase of the cycle, fluid 
enters through the outlets: in that case, velocities are put to zero (but not the mass flux) and Dirichlet 
boundary conditions are imposed for scalars (the same values as for the initialisation are used for 
turbulence quantities). 

• The domain is initialized with a low turbulence energy (1.5 x 1CT 4 m 2 s -2 ) and a dissipation rate 
giving a ratio ly /u = 10. 

Sensitivity studies 

k-e computations were performed using two different structured grids. The coarse grid is shown in Fig. 1: 

• minimum cell size (squared cells at the slot exit) : 0.0254mm 
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• maximum cell size : 12.5mm 


• number of cells : 15707 

For the fine grid, each cell has been divided into four cells, leading to 62828 cells. Some characteristics of 
the flow are sensitive to the grid, such as the width of the jet (long-time averaged flow), but the dynamics is 
by no means influenced. 

As shown in the table above, different time steps were used in order to evaluate the numerical error due 
to the time-marching scheme. The results are only marginally influenced by the reduction of the time step. 
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Solution Methodology 

Computations were carried out using a modified version of the NEAT finite volume code [1]. The convection 
terms are discretized using the essentially second order central difference based CONDIF [2] scheme. For 
the diffusion terms, central differences are used. To integrate the equations through time, the Crank-Nicolson 
scheme is used. Compatible with our future intention to model the jet with LES staggered grids are used. 

To provide a pressure field the SIMPLE algorithm [3] was used. For Unsteady Reynolds Averaged 
Navier-Stokes (URANS) solutions, the use of various solution adaptive time step methods is also being 
explored. 


Model Description 

Three URANS turbulence models were used in these simulations, low Reynolds number k-e [4], non-linear 
k-e (cubic) [5] and Explicit Algebraic Stress Model (EASM) [6]. The implementation of these models is 
described in [7]. 


Implementation and Case Specific Details 

An unsteady boundary condition was applied at the jet inlet and varied with time as shown in Eqn (1) below. 


U{t) — Asinflr ft) (1) 

Where U is the instantaneous jet amplitude (ms -1 ), A is the maximum jet amplitude of 30ms -1 , / is the 
frequency of 444.7Hz and t is the present time(sec). 

The velocity profile of this jet initially set as a top-hat profile but it was seen that applying a seventh 
power profile improved the accuracy by more correctly modelling the shear layers induced by the walls. 
Since the fluid density is assumed constant, an unsteady boundary condition was also set across the entire 
top face of the domain. To improve numerical stability, the side walls were set as solid boundaries. This is 
shown in Figure 1 

An unsteady k and e boundary condition was also set at the inlet. The turbulence intensity was assumed 
to be 10% (and isotropic) and the length scale, L , was taken as the width of the jet (1.27mm). However, on 
the suction stroke of the jet cycle, k and e were set to zero in order to simulate no turbulent energy being 
input into the domain. This is shown by Eqn (2) and (3) 


3 

k = max[-(0.1U(t)) 2 , 0] 

, P / 2 

e = max[C 3 ^ 4 — - — , 0] 

where C=0.09 and 1=0.07 L 


( 2 ) 

( 3 ) 
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Figure 1: The boundary conditions. 


The timestep was chosen to correspond to a tenth of a single degree of the cycle corresponding to 
6.25 xlO 7 s. This allowed the phase average to be taken at a discrete number of timesteps rather than 
having to resort to interpolation between timesteps. 

The grid used was of non-uniform rectangular type with 99 grid nodes in the x direction across the jet 
and 49 nodes in the streamwise y direction. The geometric grid expansion factor was typically about 1.15. 
The grid is shown in Figure 2. Although quite coarse, this grid allowed reasonably quick simulations to be 
run. However a finer grid is to be investigated in the future. 

Four iterations per timestep were found necessary for convergence. The maximum Courant number is 
1.1 (that is where the grid spacing is smallest and velocity highest) but the average Courant number is well 
below unity. Initial non-cyclic data is discarded. Typically, the flow reached a cyclic state within about 
15-20° 

The simulations were run in sequential mode on a SUNFIRE 6800 machine with 24 individual 950MHz 
processors and 12GB of working memory. This resulted in run times of roughly 60 seconds/timestep for the 
standard k-e model and 75 seconds/timestep for the cubic and EASM models. 
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Figure 2: The grid structure used. 
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Introduction 

Active flow control technologies is growing area of aerodynamic research in the early 21 st 
century. The goal is to prevent boundary layer separation and as such it is often applied to designs 
of high-lift airfoils. Forced oscillations superimposed on a mean flow which is on the verge of 
separation point is an effective means to delay boundary-layer separation, such as blowing or 
suction techniques. However, the progress in active flow control technologies has often been 
paced by the development of actuator capabilities. A popular current actuator is the synthetic jet, 
which has demonstrated capabilities regarding separation control and thrust vectoring. 

Our approach to this problem has been to employ two different computational fluid 
dynamics (CFD) codes — a two-dimensional structured grid solver and a three-dimensional 
unstructured grid solver. Both codes are in-house codes in use at the University of Kentucky. The 
structured code, GHOST, uses overlapping Chimera-style grids to model more complicated flows 
and geometries — the advantages and disadvantages of this type of approach are well-established, 
although with proper diligence excellent results can be obtained using these methods. However, 
the continual advances in both computational fluid dynamics (CFD) algorithms and computer 
technologies have made unstructured grids more attractive with their ability to smoothly conform 
to varying flow conditions and complicated boundaries with a single grid. Remaining challenges 
for this type of approach including grid generation for a computational domain with complex 
geometries, well-balanced grid decomposition on distributed system, and efficient parallel 
performance. For the three-dimensional simulation of this case, we use a three-dimensional 
unstructured grid code called UNCLE. While one is structured and one is unstructured, GHOST 
and UNCLE share some common numerical techniques and as such are quite complementary 


Solution Methodology 

The two-dimensional simulations presented were computed with GHOST. GHOST is an in-house 
CFD code developed at University of Kentucky by P. G. Huang. The code is based on a finite 
volume structured formulation with chimera overset grids. The QUICK and TVD schemes are 
applied to discretize the convective terms in the momentum and turbulence equations, 
respectively; the central difference scheme is used for the diffusive terms and the second order 
upwind time discretization is employed for the temporal terms. This code has been tested 
extensively and is routinely used for turbulence model validation [1,2,3] and flow control studies 
[4,5,6]. The turbulence model used in the present computation is Menter's SST two equation 
model [7], which provides excellent predictive capability for flows with separation [8]. The 
multi -block and chimera features of the code allow the use of fine gird patches near the jet 
entrance and in regions of highly active flow. The code also employs MPI parallelization to allow 
different computational zones to be solved on different processors. 

The unstructured grid code is UNCLE. UNCLE is a 2D/3D finite volume unstructured 
unsteady incompressible Navier-Stokes solver. UNCLE employs a pressure-based SIMPLE 


1.7.1 



algorithm with second order accuracy in both time and space. A second order upwind scheme is 
used for computing advection terms. Non-staggered grids with the Rhie and Chow momentum 
interpolation method [9] are employed to avoid checkerboard solutions. In order to take care of 
turbulence flow in most realistic cases, F. R. Menter’s shear-stress transport (SST) turbulence 
model [7] is implemented. It is designed to study the challenges of using unstructured grid codes 
on high-performance parallel computers to simulate realistic cases. To increase flexibility in 
complex geometries, UNCLE may use a variety of grid types, such as triangles, quadrilateral, 
tetrahedron and hexahedra meshes. In order to achieve good load balance for parallel computing, 
METIS [10] is used to partition the grid. Generally, there are two different partitioning 
approaches - vertex based and element based partitioning for mesh-partitioning as shown in Fig. 
1(a) and (b) respectively. For vertex-based partitioning, the boundary elements are doubled and 
the vertices at the boundary are overlapped. Since the control volumes at boundary are not 
partitioned, only communication of boundary nodal properties is required. For element-based 
partitioning, the boundary vertices are doubled. Because the control volumes at the boundary are 
split, all nodal points surround a boundary vertex are needed to interpolate the properties of the 
boundary vertices. Communication of boundary nodal properties in element based partitioning is 
heavier than vertex-based partitioning. On the other hand, vertex-based partitioning has to handle 
doubled elements at the boundary, it still costs computational time. For the purpose of direct use 
the information from METIS, element based partitioning is used for pre-processing code of 
UNCLE. 

Both UNCLE and GHOST are designed to run on multiple platforms, notably in this 
commodity PC clusters on which the bulk of the simulations were performed. 


Model Description 

The turbulence model used in the present computation is Menter' s SST two equation model!], 
which provides excellent predictive capability for flows with separation [8], among other 
challenging flow types. UNCLE likewise employs the M-SST two-equation turbulence model, 
but has not undergone the extensive “burn-in” that GHOST has received. 


Implementation and Case Specific Details 

For the current 2-D computations, we use the 2-D grid (structure 2D grid #1 and structure 2D grid 
#2) provided by the workshop website, in combination with two timesteps (T 1=0.003 and 
T2=0.009), generating 4 sets of results. The following convention is used to name the data files: 

G1T1: structure 2D grid #1, and time step Tl=0.003. 

G1T2: structure 2D grid #1, and time step Tl=0.009. 

G2T1: structure 2D grid #2, and time step T 1=0.003. 

G2T2: structure 2D grid #2, and time step Tl=0.009. 

The non-dimensional time period for one cycle is: T =52.488. 

The value being chosen to non-dimensionlize the velocity is: V = 30 ml s . 

The Reynolds number is Re=2453.72. 

The left hand side boundary of zone 1,2, 3, 4 (workshop provided grid, moving diaphragm) are 
given by oscillatory non-dimensional velocity boundary condition. The oscillation is defined by 
the displacement history at the center of the diaphragm from the PIV experiment provided on the 
cfdval2004 website translate into the velocity history at the center of the diaphragm. The time- 
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dependent velocity at the center of the diaphragm can be approximated by curve-fitting method, 
and its mathematical formulation is described as: 

v 

— = a + b- cos (c • t + d) , for (2) 

a = - 0.0044205499 
b = 0.026689782 
c = 0.119707081 
d = - 4.0776928 

where a, b, c, and d are dimensionless parameters and t is dimensionless time. The u-t plot is 
shown in Fig. 3. For the top of box, the outflow boundary condition is chosen to satisfy the 
continuity equation. For the rest of the boundary, the no-slip condition is imposed as a wall 
boundary condition. 

In the unstructured simulations, the slot width, 0.05 inches, is chosen to be the reference 
length. Air density at sea level set as 1.185kg/m 3 , viscosity is 18.4e-6kg/m-s, and the reference 
velocity is the maximum velocity at the slot exit which is approximately 30m/s. As with the 
GHOST simulations, the Reynolds number for this case is 2453.7. The single timestep used in 
these simulations is 0.02 seconds, so the timesteps per cycle is 2624. 

An unstructured grid is generated to fulfill the grid format of UNCLE by using GAMBIT. 
The total grid points are 0.26 million points and the cell are 1.4 million. For this grid, y + 
approximately equals 5 on the wall surface. The geometry of the grid is exactly the same as 
provided by the cfdval2004 website. Figure 2 (a) and (b) show the global outline and the zoom-in 
picture of the actuator, respectively. 

As with the 2D simulations, a periodic velocity boundary condition replaced the moving 
boundary condition for the diaphragm. The non-dimensional period T can be calculated by: 

1 /T = f*=f*£/u (3) 

where / is non-dimensional frequency, frequency / =444. 7Hz, reference length £ =0.05 in, and 
reference velocity u =30 m/s. As a result of equation 3, the non-dimensional period T is again 
52.48. The diaphragm boundary is approached in the same fashion as the 2D simulation (Eqn 2), 
as are the other boundary conditions. 
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Fig. l(a, left) Schematic of vertex based partitioning, (b, right) Schematic of 
element based partitioning. 
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Fig 2 (a, left) a global view of the generated unstructured mesh (b, right) 
zoom-in picture of the actuator in unstructured grid mesh 



Fig. 3 The velocity profile of periodic boundary condition in one cycle. 
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Introduction 

A synthetic jet is generated by a piezoelectric diaphragm embedded on one side of a cavity. Flow 
enters and exits the cavity through an orifice in a periodic manner. Synthetic jets have been 
shown to exert significant control authority in many applications and have the additional benefit 
of being compact with zero net mass flux. During the past several years, many experimental and 
simulation studies have been reported in the literature to characterize the behavior of synthetic 
jets and their applications in flow control. However, there has not been a systematic study to 
assess the current capabilities of different classes of turbulent flow solution methodologies to 
predict flow fields induced by the synthetic jets. The purpose of this NASA Langley workshop 
[1] is to provide a systematic evaluation of various turbulent flow simulations methodologies by 
examining their performance against three experimental test cases. 

Out of the three benchmark cases, Case 1: a synthetic jet issuing into the quiescent air, has been 
studied here by employing a 2D grid (4 zones, 35,986 grid points total). The simulations are 
performed by employing the computational code WIND- v. 5 [2], which is a multi-zone structured- 
grid compressible Reynolds- Averaged Navier-Stokes (RANS) solver. Several turbulence models 
have been tested, such as Shear-Stress Transport (SST) [5], Spalart-Allmaras (SA) [4], SST 
combined with LES. The grid size and time-step independence have also been assessed. 
Comparisons with the experimental data show that the SST model gives the best results. 

Solution Methodology 

CFD code WIND is a product of the NPARC Alliance [2,3], a partnership between the NASA 
Glenn Research Center (GRC) and the Arnold Engineering Development Center (AEDC) 
dedicated to the establishment of a national, applications-oriented flow simulation capability. 
WIND computes the solution of the Euler and Navier-Stokes equations, along with supporting 
equation sets for turbulent and chemically reacting flows by employing a variety of turbulence 
models and chemistry models respectively. WIND is coded in Fortran 77, Fortran 90, and C 
programming languages. The governing equations are solved in conservation form. Explicit 
viscous terms are computed using either upwind or central differencing, and their order may be 
controlled through the use of keywords in the input data file. In all the simulation results 
presented here, the order of accuracy is second-order. The implicit convection terms are 
computed using either an approximate factorization scheme with four-stage Runge-Kutta time- 
stepping or they may be disabled altogether, and a global Newton iteration scheme may be 
employed. WIND uses externally generated computational grids. The solution is executed 
iteratively on this grid. 
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Model Description 

As stated before, WIND is an Unsteady Reynolds-Averaged Navier-stokes (URANS) flow solver. 
WIND has various turbulence models: Spalart-Allmaras (SA) one-equation model [4] and 
Menter’s Shear Stress Transport (SST) two-equation models [5, 6]. For unsteady flows, a 
combined SST and Large Eddy Simulation (LES) [7] model is also available. The combined 
model reduces to the standard SST model in high mean shear regions (e.g., near viscous walls), 
where the grid is refined and has a large aspect ratio unsuitable for LES models. As the grid is 
traversed away from high mean shear regions, it typically becomes coarser and more isotropic; 
the combined model smoothly transitions to an LES model. The intent is to improve predictions 
of complex flows in a real-world engineering environment, by allowing the use of LES methods 
with grids typical of those used with traditional Reynolds-Averaged Navier-Stokes solvers. An 
input parameter C B (the default value 10.0 is used here) specifies the size of the RANS and LES 
zones. Increasing C B increases the size of the region in which the combined model reduces to the 
standard SST model. All these three turbulence models (SA, SST, SST/LES) have been tested in 
our simulations presented in this paper. 

Implementation and Case Specific Details 

A two-dimensional grid is generated by a pre-processing Grid MANipulation (GMAN) code [8]. 
Zonal connectivity information is also computed by GMAN, and is stored in the grid file used by 
WIND. During the course of a solution, WIND maintains continuity in flow properties across 
zone boundaries through a process known as zone coupling [9]. The grid used in our simulations 
is shown in Fig. 1. It has 4 zones and 35,986 grid points. Zone 1 (33*86) is the left side of cavity, 
zone 2 (62*50) is the right side of cavity, zone 3 (41*65) is the slot and zone 4 (197*139) is the 
external flow domain. In the external flow domain, there are 16 grid points from the wall (y=0) to 
a distance of y=0.1mm, which makes the resolution near the wall high enough to capture the 
features of the flow. Across the slot of the synthetic jet (in the x-direction), there are 41 points, 
which are enough to provide the accurate velocity profiles (w, v) across the slot. 



Figure 1: 2D grid used in the simulation 

(a) Whole view (b) Configuration of four zones (c) Zoomed in view of the grid in the slot 

All the boundary conditions are also specified in GMAN, and are stored in the grid file used by 
WIND. In the external flow domain (zone 4), the bottom wall is no-slip, from which the synthetic 


1 . 8.2 


jet ejects into the quiescent air. The other three boundaries of zone 4 (left, right and top) are 
specified as outflows. All the boundaries belonging to the cavity (zone 1, 2 and 3) are no-slip 
viscous walls, except the one where the diaphragm is located (1=1, zone 1). The “arbitrary 
inflow” is specified at this boundary as follows: 


v(x,y = const, t) = 0 

(i) 

u(x , y = const , t) = U sin co t 

(2) 

du dp 

P dt dx 

(3) 


It is clear that the “mode shape” of the piezoelectric diaphragm is not simulated. Instead, it is 
modeled as a “piston”, with velocity uniformly distributed in space. The volume change of the 
cavity is not included directly in the simulation either. It is found by numerical simulations that 
the net-mass-flux of the synthetic jet depends only on the pressure at the diaphragm, and not on 
the excitation velocity. Therefore, this pressure is tuned by calculating the mass-flux of the 
synthetic jet such that the net mass-flux is approximately zero. As shown in Fig. 2, when the 
pressure at the diaphragm is 14.709 psi, the relative net-mass-flux at the diaphragm is less than 
1.0e-4, and the mass flux at the synthetic jet slot is also less than 1%. 




Figure 2: Mass-flux at the diaphragm & SJ slot Figure 3: Pressure inside the cavity 

After this pressure is determined, the velocity amplitude of the jet at the exit depends only on the 
velocity amplitude at the diaphragm. The excitation velocity at the diaphragm is determined to be 
3.5 m/s. For incompressible flow, in Eq. (2), U=Uj b/W= 0.8m/s, where Uj is the velocity 
amplitude of the synthetic jet, and W and b are the height of the diaphragm and the width of the 
jet-orifice respectively. Fig. 2 shows that the amplitude of mass flux at the slot is less than one- 
fourth of the amplitude of mass-flux at the diaphragm. This implies that the amplitude of velocity 
at the diaphragm is about four times 0.8 m/s, i.e. 3.2m/s. This 2D simulation assumes that the 
cavity is a cube. However in the experiment, the cavity is actually a cylinder. Thus, the volume 
ratio of the cavity to the orifice is not the same as in the experimental setup. Therefore, it is 
expected that the velocity used in the simulation here is not same as in the experiment. 

Fig. 3 shows the pressure at three locations: the center of diaphragm, the center of the opposite 
side of diaphragm inside the cavity and at a location just above the jet-slot (x, y)=(0, 0.1mm). In 
the experiment, the pressure inside the cavity varies between 14.46 & 14.89 psi; it is between 
14.47 & 14.95 psi in the simulation. Fig. 4 shows the calculated v- velocity at (x, y)=(0, 0.1mm), 
which is in good agreement with the PIV & Hotwire data. These observations justify that the 
boundary conditions employed at the diaphragm in our 2D simulation are applicable to define the 
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experimental set-up, since the computed flow inside the cavity and at the slot behaves in a 
manner very similar to that in the experiment. 



Figure 4: Phase-averaged v- velocity 



Figure 5: Long-time averaged v- velocity 


at (x, y) = (0, 0.1 mm) 


at (x, y) = (0, 1 mm) 




(a) Long-time averaged v- velocity (b) Phase-averaged v-velocity at (0, 2mm) 

along the centerline 

Figure 6: Time-step and grid independence studies; dt = time-step 




(a) v-velocity along the centerline 

Figure 7: Long-time averaged v-velocity along the centerline of SJ. 
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The results presented in Figs. 2-4 have been obtained by using a time step corresponding to 
10000 steps per SJ cycle. The time-step therefore is 1/444.7 (Hz)/10000 = 2.3486e-7(sec). Time 
step & grid independence studies were conducted by using a half time step and half grid spacing 
in y direction. Fig. 6 shows the long-time averaged v-velocity along the centerline and phase- 
averaged v-velocity at (x= 0, y= 2 mm). Both Figs. 6(a) and 6(b) show the effect of time-step and 
grid refinement on the solution using a SST turbulence model. 

Different turbulence models have also been tested. Fig. 7 shows the variation of v-velocity along 
the centerline using the SST, SST_LES, and SA models, and comparison of calculations with the 
PIV and hotwire (HW) data. There is noticeable discrepancy between the PIV and HW data. Near 
the jet exit, all the simulation results lie between the PIV and HW measurements. In the range 
y= 10 to 20mm, simulated v-velocity is very close to the HW data, except for the SST combined 
with LES model, which over-estimates the velocity. The maximum velocity in the experiment is 
7.8 m/s. It is about 9 m/s in most of the simulations. The maximum velocity appears near the wall 
in the experiment, but not in simulations. Fig. 5 shows the computed long-time averaged 
v-velocity at (x, y) = (0, 1 mm). Again, there is a good agreement with PIV data. 

Based on the above comparisons, we conclude that the URANS simulation with SST model, with 
10000 time steps per synthetic jet cycle, gives the best results for this case. 
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Introduction 

The unsteady evolution of three-dimensional synthetic jet into quiescent air is studied by 
time-accurate numerical simulations using a second-order accurate mixed explicit- 
implicit fractional step scheme on Cartesian grids. Both two-dimensional and three- 
dimensional calculations of synthetic jet are carried out at a Reynolds number (based on 
average velocity during the discharge phase of the cycle Vj, and jet width d) of 750 and 
Stokes number of 17.02. The results obtained are assessed against PIV and hotwire 
measurements provided for the NASA LaRC workshop on CFD validation of synthetic 
jets. 

Numerical Methodology 

The evolution of zero-net mass-flux synthetic jet from a cavity into quiescent air is 
modeled by the unsteady, incompressible Navier-Stokes equations, written in tensor form 
as 

dii _ dii d Ujii dp 1 d 2 ii 

— = 0 , — + — — — + — 

dx ; dt dxj dx l Re dxdx, 

where the indices, i =1, 2, 3, represent the x, y and z directions, respectively; while the 
velocity components are denoted by u («i), v (uf) and w (i/3), respectively. The equations 
are nondimensionalized with the appropriate length and velocity scales where Re 
represents the Reynolds number. The Navier-Stokes equations are discretized using a 
cell-centered, collocated (non- staggered) arrangement of the primitive variables ( u,p ). In 

addition to the cell-center velocities ( u ), the face-center velocities, U , are also computed. 
Similar to a fully staggered arrangement, only the component normal to the cell-face is 
calculated and stored. The face-center velocity is used for computing the volume flux 
from each cell. The advantage of separately computing the face-center velocities has been 
discussed in the context of the current method in Ye et al. [1]. The equations are 
integrated in time using the fractional step method. In the first step, the momentum 
equations without the pressure gradient terms are first advanced in time. In the second 
step, the pressure field is computed by solving a Poisson equation. A second-order 
Adams-Bashforth scheme is employed for the convective terms while the diffusion terms 
are discretized using an implicit Crank-Nicolson scheme which eliminates the viscous 
stability constraint. The pressure Poisson equation is solved with a Krylov-based 
approach. 
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The solver uses a multi-dimensional ghost-cell methodology to incorporate the effect 
of the immersed boundary on the flow; however, the absence of curvilinear topology in 
the geometry makes this approach unnecessary in the current study. Care has been taken 
to ensure that the discretized equations satisfy local and global mass conservation 
constraints as well as pressure-velocity compatibility relations. The code has been 
rigorously validated by comparisons against established experimental and computational 
data. Details have been presented elsewhere [2, 3]. 

Implementation and Case Specific Details 

The present study models the flow inside the cavity using a pulsatile velocity boundary 
condition, v = v a sin(rur) prescribed at the bottom of the cavity (see Fig. 1) in order to 

generate a natural flow at the slot exit. The shape of the cavity is approximated to be a 
rectangular box without taking into consideration the finer details that make up the 
interior of the cavity. However, the geometrical and flow parameters used in the current 
study are chosen based on a scaling analysis of various parameters involved like slot 
width (d), slot height (h), cavity width (W), cavity height (H), the diaphragm vibrating 
frequency (f), etc.. For instance, the slot size is chosen such that the ratio h/d and W/d 
match those used in the experiments; however, H/d of 4.95 used in the calculations is not 
matched in the experiments. The Reynolds number (Re) in this work is defined based on 
average jet velocity during the discharge phase of the cycle (Vj) and jet width 
(d)i.e., Re = Vjd I v. Average jet velocity (Vj) in the numerical calculations is set equal to 

the value of 10.5 m/s, obtained by averaging the velocity provided by the LDV 
measurements during the discharge phase of the cycle. Because of the uncertainty in Vj 
reported by different measurements, a nominal Re of 750 corresponding to the lower 
bound of Vj calculated from the measurements is chosen in the present investigation. 
Experimental values of co = 2794 rad/s, d = 1.27 mm and v = 1.5527xl0" 5 nf/s are 

matched to give a Stokes number (S - cod 2 1 v ) of 17.02 in the computations. Various 
cases considered in the present study and the corresponding flow parameters are detailed 
in Table 1. 


# 

Re 

2D/3D 

S 

Exterior 
Domain Size 

Grid Size 

Time steps/cycle 
N 

, In 

At = 

coN 

1. 

750 

3D 

17.02 

30dx 30d x 3d 

132x250x 16 

14,000 

0.00116144 

2. 

750 

2D 

17.02 

30d x 30d 

132 x 220 

14,000 

0.00116144 


Table 1. Various cases considered in the study and their flow parameters. 


Figure 1 shows the schematic of the computational domain and the boundary conditions 
used in the computations. An outflow velocity boundary condition is prescribed at the 
left, right and top boundaries that allow them to respond freely to the flow created by the 
jet. In 3D calculations, periodic boundary conditions are prescribed in the span wise (z) 
direction. Figure 2 shows an x-y slice of a typical 3D mesh used in the region near the 
slot in the computations. Grids used in the current work are non-uniform in both x- and y- 
directions, and uniform in the span wise (z) direction in the case of 3D calculations. 
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Sufficient clustering is provided in the slot-region along x- and y-directions to resolve the 
vortex structures that form at the slot exit, as well as the shear layer in the slot. Typically 
32 grid points clustered using a cosine-hyperbolic distribution are used across the slot. In 
3D calculations, the three-dimensionality in the solution is instigated by introducing a 
small sinusoidal spatial perturbation in the z-component of velocity (w) over a few 
hundred time-steps in the first cycle. For all cases presented here, the first two cycles of 
calculations are not included in the averaging process to eliminate transient effects, and 
the next two to three cycles are used in the process. 

Results 

Figure 3 shows the isosurfaces of vorticity magnitude obtained from three-dimensional 
calculation before the onset of full three-dimensionality in the solution. It is clear from 
the figure that the flow is dominated by counter-rotating vortex pairs. Figure 4 shows the 
isosurfaces of vorticity magnitude for the same solution after the onset of full three- 
dimensionality. Figure 4 shows two pairs of rib vortices along the span obtained from the 
3D simulation. The plot of phase-averaged v-velocity component vs. phase angle (j) at the 
point (x, y) = (0 mm, 0.1 mm) shown in Figure 5 describes the procedure involved in 
aligning the CFD data with PIV data. Phase angle at which V avg = (V max + V min )/2 intersects 
the curve is made 340° by applying the required phase shift. Figure 6 shows the plot of 
phase-averaged v-velocity component vs. phase angle (j) at the point (x, y ) = (0 mm, 0.1 
mm) for the two calculations and it can be seen that the CFD data is reasonably aligned 
with the PIY data. Plot of time-averaged u- and v-velocity components at the point (x, y) 
= ( 0 mm, 0.1 mm) across the slot region is shown in Figure 7. Figure 8 and 9 show the 
comparison of phase- averaged u- and v-velocity components at the same point at (j) = 90° 
and (j) = 270° respectively. Three-dimensional calculations at Re > 1000 corresponding 
to the upper bound of V) reported in the experiments are being carried out and since the 
solutions at these Reynolds numbers were not converged at the time of preparation of this 
report, they would be presented in detail at the workshop. 
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Figure 1 : Schematic of the computational 
domain and boundary conditions. 



Figure 2: A typical grid in the slot region 
used in the computations. 



Figure 3: Isosurfaces of vorticity magnitude Figure 4: Isosurfaces of vorticity magnitude 
2 nd cycle, Re = 750, S = 17.02. 3 rd cycle, Re = 750, S = 17.02 



Figure 5: Plot of phase- averaged v vs. phase Figure 6: Plot of phase-averaged v vs. phase 

angle before aligning CFD data with PIV data. angle after aligning CFD data with PIV data. 
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Figure 7: Plot of time-averaged u- and v- velocities along the horizontal line 
y = 0.7874 (~lmm in the experiments). 



Figure 8: Plot of phase- averaged u- and v-velocities along the horizontal line 
y = 0.7874 (~lmm in the experiments) at </) = 90° degrees. 



Figure 9: Plot of phase- averaged u- and v-velocities along the horizontal line 
y = 0.7874 (-1mm in the experiments) at <f>= 270° degrees. 
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Introduction 

Accurate details of the general performance of fluid actuators is desirable over a range of flow conditions, 
within some predetermined error tolerance. Designers typically model actuators with different levels of 
fidelity depending on the acceptable level of error in each circumstance. Crude properties of the actuator 
(e.g., peak mass rate and frequency) may be sufficient for some designs, while detailed information is needed 
for other applications (e.g., multiple actuator interactions). 

This work attempts to address two primary objectives. The first objective is to develop a systematic 
methodology for approximating realistic 3-D fluid actuators, using quasi-l-D reduced-order models. Near 
full fidelity can be achieved with this approach at a fraction of the cost of full simulation and only a mod- 
est increase in cost relative to most actuator models used today. The second objective, which is a direct 
consequence of the first, is to determine the approximate magnitude of errors committed by actuator model 
approximations of various fidelities. This objective attempts to identify which model (ranging from simple 
orifice exit boundary conditions to full numerical simulations of the actuator) is appropriate for a given error 
tolerance. 


Solution Methodology 

The solution methodology used for this work is described in detail elsewhere [1]. Only the general aspects 
of the approach will be replicated herein. The time-dependent 2-D Navier-Stokes equations are used to 
describe the unsteady compressible flow generated by a synthetic jet actuator. No turbulence model is used 
in these simulations. The governing equations in curvilinear coordinates (£, rf) are written in conservation 
law form. All the length scales and dependent variables have been nondimensionalized by the orifice width 
d and the corresponding reference values, respectively, except for p which has been normalized by p oo^, 
where has been chosen to be one tenth of the speed of sound. The viscosity coefficient p is assumed to 
be constant, and the equation of state for a perfect gas is used to relate pressure to the conservation variables. 

The governing equations are closed with the following boundary conditions. A no-slip boundary condi- 
tion for the velocity vector and a constant wall temperature are imposed on the wall surface 

U \wall = V \wall = 0 ’ T \wall = ( 1 ) 

At the subsonic outflow boundary, a boundary condition for the pressure is imposed weakly. Characteristic 
conditions are applied at the upper boundary so that the vortex structures can leave the computational domain 
without producing perceptible spurious reflections. The unsteady flow inside the actuator cavity, generated 
by harmonic motion of the diaphragm, is modeled by using a new reduced-order model described in the next 
section. 

The actuator geometry is represented as the summation of multiple subdomains, each having a high de- 
gree of smoothness. Within each subdomain, a computational grid is generated having sufficient smoothness 
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Figure 1: The 2-D and quasi- 1-D geometries. 


to allow solutions with three significant digits of accuracy. The connectivity at the subdomain interfaces is 
Co smooth: thus meshlines connect at interfaces but need not have smooth derivatives. 

Two distinct geometries are the focus of the study. The first is a simple 2-D projection of the 3-D 
actuator used in Case I. The 2-D geometry of the actuator is not the centerline cut of the actual geometry, 
but rather is adjusted slightly so that the volume of the 2-D actuator (assuming a constant width of 1.24”) 
is precisely that of the 3-D actuator. This was done because the volume of the actuator strongly affects 
the resonance characteristics of the device. The exterior portion is truncated at 100 orifice diameters. The 
second geometry is a subset of the first, but only includes the exterior portion of the device and the nozzle 
portion of the actuator. The rest of the actuator is approximated with the reduced order model, by solving 
the quasi- 1-D Euler equations. The volume of the quasi- 1-D actuator has also been chosen to be equal to the 
corresponding volume of the 3-D actuator. Both the 2-D and quasi- 1-D geometries are shown in Fig. 1. The 
region where the quasi- 1-D model is used is bounded by the dashed line, whereas in the rest of the domain, 
the 2-D Navier-Stokes equations are solved. 

The grid was adjusted within the actuator and in the nearfield of the exterior to achieve two significant 
digits in the solution accuracy. The exterior grid was sufficient in resolution to accurately resolve the vortices 
near the orifice and up to ten diameters away from the orifice. Beyond that point, the mesh was expanded 
and the vortices were allowed to diffuse. The total number of grid points in the full 2-D and quasi- 1-D 
simulations were 60,697 and 38,805, respectively. 

A fourth-order upwind-biased linear finite difference scheme based on the local Lax-Friedrichs flux 
splitting is used to discretize both the 2-D Navier-Stokes equations and the quasi- 1-D Euler equations. For 
sufficient grid smoothness within each 2-D subdomain, design order accuracy is achieved. The interfaces 
are connected using the penalty approach (SAT) described in detail in references [2, 3, 4, 1]. The interface 
penalty treatment is conservative, and maintains the underlying accuracy of the interior scheme. The SAT 
procedure requires co-located solution variables on both sides of the interface. The difference between the 
two solution values was used as a measure of spatial grid resolution. For all the simulations, the interface 
error was less than 0.5%. 

The semi-discrete equations were explicitly integrated in time with a low-storage 4th-order Runge-Kutta 
scheme [5]. The simulations were all run at the maximum stable timestep. The temporal scheme has an error 
estimator that monitors the temporal error per timestep. The timestep error varied over the period in the range 
lO" 10 — 10 -7 , based on the L ^ norm of the density variable. Other solution variables had similar error 
norms. Each cycle required about 10 5 timesteps to complete. 


1 . 10.2 



Model Description 


A gap exists between the 0-D models and the full 2-D/3-D models of a synthetic jet actuator. To combine 
the accuracy and conservation properties of the full numerical simulation methods with the efficiency of 
the simplified blowing/suction type boundary conditions, a new reduced-order model of a multidimensional 
synthetic jet actuator is proposed. In contrast to the methods available in the literature, the new approach uses 
a reduced-order model to approximate a 2-D or 3-D actuator. The multidimensional actuator is simulated 
by solving the time-dependent quasi- 1-D Euler equations. The time-dependent quasi- 1-D Euler equations 
can be written in the following conservation law form: 
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where A is the cross-sectional area of the quasi- 1-D actuator. It is assumed that A is a continuously differ- 
entiable function that is independent of time, i.e., A = A(y). 

To simulate the diaphragm dynamics, a time-dependent one-to-one coordinate transformation, 


C = «*,»), <3) 

is employed to map a physical domain with the moving boundary onto a unit interval. Note that the ( 
coordinate depends on time and, therefore, a moving mesh technique is applied to solve the quasi- 1-D Euler 
equations. Because the frequency of diaphragm oscillations a; is a given quantity, the moving mesh can be 
generated analytically 

y((, t) = (1- C) [L + a( 1 - cos (lot ))] , (4) 

where y and (' are physical and computational coordinates, respectively, a and u are the amplitude and 
frequency of diaphragm oscillation, and L + a is a mean depth of the quasi- 1-D synthetic jet actuator. 

Diaphragm oscillations are forced by varying the position of the diaphragm y( 0, r) where the imper- 
meable wall boundary condition is imposed. Because the deforming mesh Eq. (4) is given analytically, the 
diaphragm velocity can be calculated by differentiating Eq. (4) with respect to time to give 

i;(0, t) = auism^Lur). (5) 

It should be noted that the region near the jet exit requires special consideration. The full numerical 
simulation of the actuator orifice region is crucial for accurate prediction of the interaction between the syn- 
thetic jet and the external boundary layer. This region is characterized by strong flow separation that cannot 
be described by the quasi- 1-D Euler equations. To overcome this problem, the quasi- 1-D actuator model 
is used only to simulate the flow inside the actuator cavity, while the small region near the actuator orifice 
is modeled by solving the 2-D unsteady Navier-Stokes equations. This approach allows us to accurately 
predict the interaction of the synthetic jet with the external boundary layer and to resolve vortices generated 
in the vicinity of the actuator orifice, while reducing the computational cost. 

The low-dimensional actuator model has several advantages. First, this approach is fully conservative 
and provides conservation of mass, momentum, and energy. Second, the new quasi- 1-D model is computa- 
tionally much more efficient compared with the 2-D or 3-D numerical simulation of the cavity flow. Third, 
the reduced-order model retains some important multidimensional features of the realistic actuator, such as 
the length, diaphragm deflection, and area variation. These properties of the new model and its ability to 
account for the compressibility effects inherent in actuator devices make it an efficient tool for quantitative 
study of the actuator resonance characteristics. 

To demonstrate the ability of the new reduced-order model to quantitatively predict the 3-D actuator 
dynamics, the phase-averaged time history of vertical velocity over the center of the slot obtained with the 
quasi- 1-D model and three different measurement techniques is shown in Fig. 2. 
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-V hotwire data, y=0.1 mm 



Figure 2: Phase-averaged time-history of the vertical velocity over the center of the slot, obtained with the 
quasi- 1-D model and using the PIV, hotwire, and LDV measurement techniques. 


Implementation and Case Specific Details 

The simulations were started with quiescent flow and run through sufficient periods to obtain periodic so- 
lutions. This was usually about 10 periods. The experimentally determined diaphragm position was used 
to describe the moving boundary condition of the quasi- 1-D actuator. The shape of the diaphragm was as- 
sumed to be the first mode of the Bessel function eigensolution of the cylindrical dram symmetric vibration. 
A scaling factor of 1.35 was later added to obtain the peak exit velocity that was in the Mach = 0.1 range. 
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Introduction 

The approach described here is to model the response of a synthetic jet using lumped element 
modeling to complement more rigorous and expensive numerical simulations. 

Solution Methodology 

The approach used is to model the actuator orifice impedance without a grazing boundary layer and is 
based on a lumped element modeling (LEM) technique, following the recent paper by Gallas et al. [1]. 

Model Description 

In LEM, the individual components of a synthetic jet are modeled as elements of an equivalent electrical 
circuit using conjugate power variables (i.e., power = generalized “flow” x generalized “effort” 
variables). Figure 1 shows an equivalent circuit representation of a piezoelectric-driven synthetic jet 
actuator, where the lumped parameters represent generalized energy storage elements (i.e., capacitors and 
inductors) and dissipative elements (i.e., resistors). Model parameter estimation techniques, assumptions, 
and limitations are discussed in Gallas et al. [1]. The frequency response function of the circuit is derived 
to obtain an expression for Q out /V ac , the volume flow rate during the expulsion part of the cycle per 
applied voltage. LEM provides a compact nonlinear analytical model and valuable physical insight into 
the dependence of the device behavior on geometry and material properties. 



Figure 1 : Equivalent circuit model of a piezoelectric-driven synthetic jet actuator. 

With slight modifications, the model is applicable to any type of zero-net mass flux actuator, 
(piezoelectric, voice-coil, etc.) and also for two-dimensional slots and axisymmetric orifices. Case 1 is 
similar to the device studied and modeled by Gallas et al. [1], namely a piezoelectric-driven synthetic jet 
exhausting into a quiescent medium through a circular orifice. 
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The various lumped parameters shown in Figure 1 are fully described in Gallas et al. [1]. The only 
empirical constant is the structural damping ratio for the acoustic diaphragm resistance, R aD . Simple 

analytical expressions are available for all other elements. The only difference for the current case from 
that considered in [1] comes from the expressions for the orifice impedance due to the use of a rectangular 
slot. Using the same approach described in [1], the acoustic resistance of the orifice is obtained for low 
frequencies assuming fully-developed laminar channel flow in a slot of width d , height h and length w 


_ V L _ 

Q out 2w(d/2f ’ 


( 1 ) 


where p is the dynamic viscosity of the fluid and Q out is the volume flow rate produced by the 
differential pressure A P out across the slot. The acoustic mass in the slot is obtained by integrating the 
distributed kinetic energy and equating it to the lumped kinetic energy, 


r _ 3 ph 
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The acoustic radiation mass M aRad is modeled for kd < 1 as a rectangular piston in an infinite baffle by 
assuming that the rectangular slot is mounted in a plate that is much larger in extent than the slot size, 


X aRad =J^ aRad 

dw 
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( 3 ) 


where X aRad corresponds to the acoustic radiation reactance, k = oo/c 0 is the wave number, and c 0 is the 
isentropic speed of sound, co is the radian frequency, and j = . 


Implementation and Case-Specific Details 

LEM provides the frequency response of the actuator for a given input voltage. Since the orifice/slot 
acoustic resistance is nonlinear, the frequency response is a function of the input voltage amplitude. As 
discussed in Gallas et al. [1], a piezoelectric-driven synthetic jet actuator exhibits two resonant 
frequencies. A definition of the various frequencies used hereafter is given below: 

/exp = 444.7 Hz = drive frequency of the diaphragm used in the experiment 

f D = 460.2 Hz = calculated natural frequency of the diaphragm, determined from linear composite 
plate theory [2] 

f H = 191 1.7 Hz = Helmholtz resonant frequency of the cavity ^f H = \j 2n^M aN + M aRad ) C aC j 

f x = 445.1 Hz = calculated first resonant peak of the entire system 
f 2 = 1976.9 Hz = calculated second resonant peak of the entire system 

The relationship between the frequencies is given by f\ f 2 = f D f H , but f\ ^ f D and f 2 ^ f H . In order to 
successfully model the piezoelectric-diaphragm that drives the actuator, the exact dimensions and 
material properties are required as shown in Figure 2. Note the large uncertainties/tolerances associated 
with the dimensions of the device. 
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Figure 2: Schematic of the Murata type 7BB-50M-1 disk used in Case 1 


Thus, in order to calculate the resonant frequency of the diaphragm we used the linear composite plate 
model developed by Prasad et al. [2]. Table 1 summarizes the input and output parameters of the 
piezoceramic driver. 


Shim (Brass) 

Elastic Modulus (Pa) 

8.963xl0 10 

Poisson’s Ratio 

0.324 

Density (kg/m 3 ) 

8700 

Thickness (mm) 

0.15 

Diameter (mm) 

49 


Piezoceramic (PZT-5A) 


Elastic Modulus (Pa) 

6.3xl0 10 

Poisson’s Ratio 

0.33 

Density (kg/m 3 ) 

7700 

Thickness (mm) 

0.31 

Diameter (mm) 

24.4 

Rel. Dielectric Constant 

1750 

d 3l (m/V) 

-1.75xlO' 10 

C EF (F) 

1.12 xlO" 7 

Voltage (V) 

101.8 


Calculated Lumped Parameters 


C aD (S 2 m4/k §) 

5.8 xlO' 11 

M aD (kg/m 4 ) 

2.062xl0 3 

t (Pa/V) 

2.2 


Table 1: Piezoelectric-diaphragm specifications and parameters. 


The measured amplitude of the phase-averaged vertical velocity component (provided in the workshop 
website) is used for comparison with the model. Three different measurement techniques were used, but 
the experimental data chosen for comparison here are the LDV measurements performed at -0.1 mm 
above the center of the slot. The experimental value is referred to as U CL (exp) , while the model output is 
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denoted as U CL . The input ac voltage used is 101.8 V, which is the mean of the fluctuating voltage given 
as input data (+100.97 V to -102.58 V). 

The lumped element model predicts (with the only assumed empirical parameter being ^ = 0.019, the 
diaphragm damping coefficient) the normalized response, U CL /U CL (exp), shown in Figure 3. In this 

case, the first resonant peak is dictated by the natural frequency of the diaphragm, and the second 
resonant peak is governed by the Helmholtz frequency of the cavity. 



Frequency [Hz] 

Figure 3: Frequency response of the synthetic jet actuator 


The results are summarized in the following table: 



V, ( ex P) Ms] 

V:l (./,) [m/s] 

Exp. (LDV) 
LEM 

24.8 

24.7 


Table 2: Comparison between experiment and model. 


Clearly, the model results are in close agreement with the experimental measurements. However, the 
important point is not whether or not the model matches the experiment, but that the model gives a 
reasonable estimate (typically within ±20%) with minimal effort. The power of LEM is its simplicity and 
its usefulness as a design tool. LEM can be used to provide a reasonable estimate of the frequency 
response of the device as a function of the signal input, device geometry, and material and fluid 
properties. 
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Introduction 

The second case for this workshop builds upon the isolated synthetic jet of Case 1 by adding a crossflow, 
with no streamwise pressure gradient, for the developing jet to interact with. Formally, Case 2 examines 
the interaction of a single, isolated, synthetic jet and a fully turbulent zero-pressure gradient boundary layer. 
The resulting flow has many of the characteristics that need to be modeled with fidelity if the results of the 
calculations are to serve as the basis for research and design with active flow control devices. These include 
the turbulence in the boundary layer, the time-evolution of the large vortical structure emanating from the 
jet orifice and its subsequent interaction with and distortion by the boundary layer turbulence, and the effect 
of the suction cycle on the boundary layer flow. 

In a synthetic jet, the flow through the orifice and out into the outer flowfield alternates between an exhaust 
and a suction cycle, driven by the contraction and expansion of a cavity internal to the actuator. In the 
present experiment, the volume changes in the internal cavity are accomplished by replacing one of the rigid 
walls of the cavity, the wall opposite the orifice exit, with a deformable wall. This flexible wall is driven by 
a bottom-mounted moveable piston. The piston is driven electro-mechanically. The synthetic jet issues into 
the external flow through a circular orifice. In the present experiment, this orifice has a diameter of 0.250 
inches (6.35 mm). The flow is conceptually similar to that documented in Schaeffler [1]. 

To document the flow, several measurement techniques were utilized. The upstream boundary conditions 
(in-flow conditions), and several key phase-averaged velocity profiles were measured with a 3-component 
laser-Doppler velocimetry system. Phase-averaged velocity field measurements were made with both stereo 
digital particle image velocimetry and 2-D digital particle image velocimetry as the primary measurement 
system. Surface pressure measurements were made utilizing an electronically scanned pressure system. 


Facility and Actuator Details 


This research effort was conducted in the NASA Langley 15-inch Low Speed Wind Tunnel of the Flow 
Physics and Control Branch. This tunnel is a closed-return atmospheric facility dedicated to basic flow 
physics research efforts. The tunnel has a maximum speed of 115 ft/s and a turbulence intensity of less 
than 0.13%. The tunnel medium is air at sea level conditions. During this research effort, the tunnel was 
equipped with a suspended flat plate model that acts a splitter plate and facilitates the installation of the 
actuator within the plate model. The plate features an elliptical leading edge. Immediately downstream of 
the transition from the leading edge to the plate, there is a grit strip to trip the boundary layer of the plate 
and also grit strips at corresponding locations on the tunnel walls. The ceiling of this tunnel is adjustable 
at several locations down the length of the tunnel. The ceiling height was adjusted to create a zero pressure 
gradient over the plate starting from a station 24.8 inches (630 mm) downstream from the leading edge to 
a station 58.8 inches (1494 mm) downstream from the leading edge. The centerline of the jet was located 
38.8 inches (986 mm) downstream of the leading edge. This results in the region of zero pressure gradient 
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Figure 1 : Pressure distribution down the centerline of the tunnel. X is measured in the streamwise direc- 
tion, and X = 0.00 is the origin of the model coordinate system. The model coordinate system is shown 
schematically in Figure 4. 


extending 56 jet diameters upstream of the jet centerline and 80 jet diameters downstream, as shown in 
Figure 1. These pressures were measured via a row of static pressure taps which are installed along the 
centerline of the tunnel with a total of 15 taps upstream of the jet exit and 15 taps downstream of the jet 
centerline, or exit. A set of span wise pressure taps, 16 in total, were located 8 jet diameters upstream of the 
jet centerline and another set of spanwise taps, also 16 in total, were located 8 jet diameters downstream. 
All of these taps were connected via short lengths of flexible Tygon tubing to two ESP scanning pressure 
transducers, each with a full scale range of ±10 inches of water. The ESP modules were mounted inside the 
plate model and the measured pressures are referenced to the static pressure of the tunnel, as observed by a 
Pitot-static tube upstream of the splitter plate. The modules were automatically calibrated at regular intervals 
during a given test day. The pressure data was used to calculate the freestream velocity over the plate and to 
verify that the pressure gradient over the plate, in the area of interest, was zero. Flow visualization conducted 
in this tunnel suggests that the flow is two-dimensional in the zero pressure gradient region of the plate, with 
the exception of the region close to the side walls [3]. 

The synthetic jet utilized in this research effort featured a circular orifice 0.25 inches (6.35 mm) in diameter. 
The throat of the orifice is smoothly tapered from a diameter of 0.60 inches (15.2 mm) on the inside cavity 
wall to the 0.25 inch diameter exit dimension. The actuator design is shown schematically in Figure 2. 
The active element of the actuator was the wall of the cavity opposite the exit of the jet, the bottom wall 
in Figure 2. This wall was displaced electro-mechanically by a sinusoidal voltage. To document the in 
situ performance of the actuator, the displacement of the actuator diaphragm was recorded, as was the 
cavity pressure and temperature. Two dynamic pressure transducers were installed in the interior of the 
cavity, one transducer referenced to the static pressure of the tunnel upstream of the splitter plate and the 
other transducer is an absolute transducer. A thermocouple was also located in the cavity. The diaphragm 
displacement and cavity pressure as a function of phase can be seen in Figure 3. The presentation of all 
of the subsequent velocity measurements utilize this same phase reference. The diaphragm displacement is 
relative to the nominal cavity depth when the tunnel is on condition. This value is 1.68 mm. 


Boundary Conditions 


In addition to the boundary conditions associated with the actuator, i.e., the diaphragm displacement, cavity 
pressure and temperature, the boundary conditions for this case additionally consist of the upstream, or in- 
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Figure 2: Schematic of the Synthetic Jet Actuator. 



Figure 3: Diaphragm displacement, drive signal and cavity pressure as a function of phase for the Synthetic 
Jet Actuator. 


flow, boundary condition, and the geometry of the model and the tunnel itself. The tunnel test condition 
was a Mach number in the test section over the plate of M = 0.10. The tunnel dimensions at the test section 
are 15.0 inches (381 mm) wide by 9.8 inch (249 mm) high (distance from the splitter plate to the top wall). 
The orifice diameter was 0.25±0.005 inches. The location of the jet orifice down the plate is specified in 
reference to the model coordinate system. The model coordinate system is shown schematically in Figure 4 
and its origin was established where the first line of spanwise pressure taps intersects the line of centerline 
pressure taps. This location is 36.8 inches (935 mm) from the leading edge of the plate. In the model 
coordinate system, the centerline of the jet orifice is located at 2±0.020 inches (50.8±0.50 mm ). Moving 
from the origin of the model coordinate system to the jet centerline in the streamwise direction defines the X 
coordinate axis, the Y coordinate axis is in the spanwise direction and the Z coordinate axis is in the vertical 
direction. The corresponding velocity components will be U, V, and W. From this point forward, distances 
in the streamwise, X, direction will be stated in millimeters or at times jet diameters from the jet centerline. 
The key locations are: x = 44.45 mm ( -ID), x = 50.80 mm ( 0D, CL), x = 57.15 mm ( ID), x = 63.50 mm ( 
2D), x = 69.80 mm ( 3D), and x = 76.20 mm ( 4D). 

While the atmospheric conditions varied on a daily basis, the conditions were never far from standard atmo- 
spheric conditions at sea level in a wind tunnel vented to the atmosphere, in a temperature-controlled room. 
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Figure 4: Schematic of the Model Coordinate System. 


Freestream Mach Number 
Atmospheric Pressure 
Atmospheric Temperature 
Density 
Viscosity 


M oo 

0.10 

P atm 

101,325 Pa 

T 

oo 

23.9 °C (75 °F) 

Poo 

1.185kg/m 3 

poo 

18AxlO~ 6 kg/ms 


Table 1 : Atmospheric Conditions Observed for Case 2 


These conditions are given in Table 1. 

The upstream boundary conditions were measured at station x = 0.000 in the model coordinate system. A 
complete description of the state of the turbulent boundary layer was measured using the 3 -component laser- 
Doppler velocimetry system to be described in the following section. All six components of the turbulent 
stress tensor were documented and presented as the upstream boundary condition. The velocity means and 
the turbulent stress components were compared to similar values available in the open literature, as shown 
for the mean profile of the freestream component of the velocity, U, in Figure 5. 


Instrumentation and Data Acquisition 

The velocity measurements to support Case 2 were carried out using a 3-component laser-Doppler velocime- 
try (LDV) system, a stereo (3D) digital particle image velocimetry (SDPIV) system, and a 2D digital particle 
image velocimetry (2D-PIV) system. The two PIV systems were the primary technique for measuring the 
velocity field. The LDV system was an orthogonal, crossed-fringe, fiber-optic probe configuration, with the 
probes mounted 90° from one another. The 514.5, 496.5, and 476.5 nanometer wavelengths from an Argon- 
Ion laser were used to measure the spanwise (V), vertical (W), and streamwise (U) velocity components, 
respectively within the tunnel. Bragg cells were utilized to provide directional measurement capability in 
all three velocity components. Both fiber optics probes used 750-mm focal length lenses, which along with 
an input beam diameter of 3 to 4mm, generated a sample volume calculated to be approximately 100 fim in 
diameter and spherical in shape. The optical setup of the LDV system can be seen in Figure 6. All velocity 
measurements were made in coincidence mode. Coincidence mode requires that a validated velocity signal 
be made independently in each of the three components within a very short time window, guaranteeing that 
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Figure 5: Upstream Boundary Condition on U compared to the data of Klebanoff[4] 


the velocity signals come from a single particle traversing the measuring volume. The use of coincidence 
mode is required to make full turbulence measurements with an LDV system. The use of coincidence in the 
present research effort allows each of the six components of the turbulent stress tensor to be measured. 

Additionally, the actuator drive and sync signal were acquired with each trigger and stored if the point was 
validated as being in coincidence. Consideration of the actuator drive and sync signal allows the LDV data, 
which is acquired at random intervals, to be reconstructed as a function of phase. The fiber optics probes 
were mounted on a X95 rail frame traversing system. The traversing system utilized several motorized slides 
that provided 1 meter of travel in each of the three axes, with 1 /i m resolution in each direction. For the 
LDV portion of the experiment, the flow was seeded with mono-disperse, 0.86 fim polystyrene latex (PSL) 
micro-spheres. The seed particles were suspended in 100-proof alcohol and atomized by a six-jet atomizer. 
The atomized particles were then introduced into the settling chamber of the tunnel. The particles were 
fabricated at NASA Langley using the technique described by Nichols [2]. 

In processing the data from the LDV system, we wish to relate the velocities measured to a particular phase 
of the actuator. The LDV system acquires, along with the coincident velocity data, data for the drive signal 
of the actuator and a sync signal. Analysis of these values allows for the reconstruction of the phase angle 
at the time of acquisition and was used to phase average the LDV data record at every point. Once the phase 
angle of each piece of LDV data is calculated, the data can be grouped into bins according to phase angle. 
Each bin represents data with a phase angle in a certain range. For the data in Case 2, this range was 10°. 
The bin was then averaged to yield the phase-average at that point. The uncertainty in the phase-averaged 
LDV mean data is estimated to be ±0.5%. Once all the phase-averaging was accomplished, each of the 
phase-averaged series was averaged as a time series to calculate the long-time average at each point. The 
long-time average velocity profiles derived from the LDV data are shown in Figure 7. 

The timing control, image acquisition, data management and post-processing of both the stereo PIV and 
the 2-D PIV data was handled by a commercial system. The light source utilized was a pulsed, frequency- 
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(a) Laser and Optical Configuration 


(b) Fiber Optic Bundles and 
Traversing System 


Figure 6: The Laser-Doppler Velocimetry Setup. 


doubled 300 mJ Nd-YAG laser operating at 10 Hz. For the PIV portion of the experiment, the flow within 
the tunnel was seeded with atomized mineral oil with a typical particle size of 2 microns. The seed was 
injected into the flow in the tunnel settling chamber. Standard sheet-forming optics were utilized and the 
cameras were positioned on opposite sides of the tunnel for the SDPIV data. Also for the SDPIV, the 
camera mounts allow for the required tilting of the lens, relative to the camera body, in order to satisfy the 
Scheimpflug condition. The light sheet orientation for the SDPIV data was in the spanwise direction and 
in the streamwise direction for the 2D-PIV data. The position of the sheet in the spanwise direction was 
controlled by a prism mounted on the same traversing system utilized by the LDV system. The light sheet 
thickness was typically 1.5 mm thick. The resolution of each camera was 1008 x 1018 pixels and the typical 
field of view was 40 mm x 40 mm for the SDPIV and 18 mm x 18 mm for the 2D-PIV. 

The image acquisition and lasers for both PIV systems were synchronized with the drive signal of the 
actuator so that images could be acquired at a specific time delay after the beginning of a cycle of the 
actuator drive signal. A specific delay for the laser pulsing and image acquisition corresponds to a specific 
phase of the actuator drive cycle. The same sync signal that the LDV system acquired was used as the 
trigger for the PIV signal. Within the PIV software, the timing delay was varied in order to acquire data at 
36 specific phases of the synthetic jet actuator drive signal, equally spaced over both the exhaust and suction 
cycles. The flowfield over the entire drive signal was recorded as 10,800 instantaneous velocity fields, 300 
velocity field measurements at each of 36 specific delays, which results in the measurement of the velocity 
field every 10 degrees of phase. The 300 image pairs were acquired over a period of 3 minutes and were 
processed using a standard PIV cross-correlation data reduction technique and then averaged to give the 
phase-averaged, or ensemble-averaged, velocity field. The uncertainty in the phase-averaged 2D-PIV mean 
data is estimated to be ± 1.3%. These phase-averaged velocity fields were be ordered in time allowing the 
evolution of the jet-induced velocity field, and its interaction with the crossflow, to be studied. The long- 
time average of the velocity field was calculated by averaging, through one complete cycle, the individual 
phase-averaged velocity fields, as shown in Figure 8 for data derived from the 2-D PIV system. 
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Figure 7: Long-Time Average Velocity Profiles from the LDV data. 




(a) Velocity Vectors (b) Velocity Vectors and Streamlines 

Figure 8: Long-time average Velocity Field via 2-D PIV. In these figures, x = 0 corresponds to 
X = 50.325 mm in the model coordinate system. 
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CASE 2: DETAILS AND SUBMISSION GUIDELINES 
Relevant Details 

The freestream Mach number is M_freestream = 0.10. The atmospheric conditions varied, 
but were essentially standard atmospheric conditions at sea level, in a wind tunnel vented to the 
atmosphere, in a temperature-controlled room. These conditions can be given as approximately: 
p_ambient = approx 101325 kg/(m-s A 2) T_ambient = approx 75 deg F (approx 297 K). Some 
derived relevant conditions are: density_ambient = approx 1.185 kg/m A 3 viscosity_ambient = 
approx 18.4e-6 kg/(m-s), u_freestream = approx 34.6 m/s, and Re_freestream = approx 2.23e6 
per meter. 

The upstream boundary conditions from the experiment (associated with the boundary layer 
on the plate at location x=50800 microns (50.8 mm) upstream of the center of the jet orifice), to 
be used to help set/verify CFD inflow BCs, are given in the file located on the website. 

The diaphragm frequency = 150.0 Hz. The "neutral position" of the moveable diaphragm 
plate at rest inside the cavity is approximately 2.8 mm below the upper cavity wall (see 
Geometry page for details). However, when the wind tunnel is on and when the moveable piston 
is operating, the neutral position of the plate moves up so that it is approximately 1.7 mm below 
the upper cavity wall. This change in the neutral position is not accounted for in the current 
CFD grids available from this website. However, it is not known whether accounting for this 
shift is important or not. The plate then displaces sinusoidally about this new neutral position 
with a maximum displacement of approximately +-0.77 mm. The file located on the website lists 
data from the cavity, including cavity pressure, voltage input, and displacement data about its 
neutral position. 

The pressure is measured inside the cavity on its top wall, to one side of the orifice, 
approximately on the order of 10 mm from the side wall. The pressure is measured with respect 
to the tunnel static pressure. The following figure shows the measured phase- averaged 
streamwise and vertical velocities over the (approximate) center of the orifice as a function of 
phase. 


Measured streamwise and vertical velocities just 



2 . 2.1 



These LV-derived data over the orifice (plus other quantities of interest) are given in a file 
located on the website. Note in the data file that the spanwise (v) velocity as measured is not near 
zero. 

Submission Guidelines 

Numerical predictions of this type of statistically unsteady flow are relatively new. The 
purpose here is to determine the state-of-the-art in modeling these types of unsteady synthetic- 
jet-type flows, so we want to explore which CFD methods work and which do not. 

• Either model the internal cavity or specify an unsteady boundary condition at the orifice 
exit. 

• Either model the moving diaphragm, or you may employ an unsteady boundary condition 
that approximates its effect within the cavity. 

There is the requirement that you detail specifically how you choose to model the case, 
including all boundary conditions and approximations made. As we assess the methodologies 
used at the workshop, it will be important to know as many details as possible about the 
calculations/ simulations. 

Detailed requirements include: 

1. The case must be run time-accurately and in three-dimensions, in order to simulate the 
unsteady 3-D nature of the case. 

• la. RANS codes run in time-accurate mode (e.g., URANS) solve directly for phase- 
averaged variables, i.e. <Uj> and <u'iUj>. (see Appendix in Case 1: Details and Submission 
GuidelinesClick). Therefore, RANS simulations should result in repeating or very-nearly- 
repeating periodic solutions. When periodicity is achieved, averages over one or more periods of 
oscillation yields the long-time-averaged (time-independent mean) values for these quantities. 

• lb. DNS, LES, or DES simulations will need to be post-processed to obtain both the 
phase-averaged and long-time-averaged (time-independent mean) values. 

2. GRID STUDY: Solutions using more than one grid size are encouraged, but not required. If 
you use more than one grid, submit each set of results separately. 

3. TIME STEP STUDY: Solutions using more than one time step are encouraged, but not 
required. If you use more than one time step, submit each set of results separately. 

Specific quantities that result from your computations at particular locations will be required 
for submission. Note that for all the following, we adopt the coordinate system with x 
downstream, z up, and y spanwise, with the (x,y,z)=(0,0,0) origin on the tunnel splitter plate 8 
diameters (50.8 mm) directly upstream of the center of the orifice (the orifice diameter is 0.25 
inches = 6.35 mm). The requirements follow (if it is not possible to provide a particular quantity, 
simply leave it out of the "variable" list, and reduce the number of columns of data submitted): 
a. Long-time-averaged downstream velocity (u), spanwise velocity (v), and vertical velocity 
(w) profiles (nondimensionalized by Uinf) along vertical lines at the centerplane (y=0) and: x=0 
mm, x=44.45 mm (-1D upstream), x=50.8 mm (center of orifice), x=57.15 mm (ID 
downstream), x=63.5 mm (2D downstream), x=76.2 mm (4D downstream), and x=101.6 mm 
(8D downstream). Give these data to at least a height of 50 mm. Also, submit horizontal lines of 
results at x=57.15 mm and: z=5 mm, z=10 mm, and z=20 mm; and also at x=63.5 mm and: z=5 
mm, z=10 mm, and z=20 mm. Give these data to at least a width of 25 mm to either side of the 
orifice. Also, submit horizontal lines of results at y=0, z=0.4 mm over the slot (from approx 
x=47.625 mm to x=53.975 mm), and at y=0, z= 1 0 mm from at least x=45 mm to x=70 mm. 
Name this file: case2.avgvel.ANYTHING.dat-where "ANYTHING" can be any descriptor you 
choose (should be different for each file if you are submitting multiple runs)-the file should be in 
6-column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #number of time steps per cycle (pound sign needed) 

6th line: #brief description of code/method (pound sign needed) 

7th line: #other info about the case, such as spatial accuracy (pound sign needed) 
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8th line: #other info about the case, such as turb model (pound sign needed) 

9th line: #other info about the case (pound sign needed) 

10th line: variables="x, mm","y, mm ","/., mm","u/Uinf","v/Uinf","w/Uinf" 

1 1th line: zone t="data along x=0, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along x=0, 
y=0, next line: zone t="data along x=44.45mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=44.45, y=0, next line: zone t="data along x=5().8mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=50.8, y=0, next line: zone t="data along x=57.15mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=57.15, y=0, next line: zone t="data along x=63.5mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=63.5, y=0, next line: zone t="data along x=76.2mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=76.2, y=0, next line: zone t="data along x=l 01 ,6mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=101.6, y=0, next line: zone t="data along x=57. 1 5mm, z=5mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=57.15, z=5, next line: zone t="data along x=57.15mm, z=10mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=57.15, z=10, next line: zone t="data along x=57.15mm, z=20mm" subsequent lines: 
x(irnn) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along x=57.15, z=20 next 
line: zone t="data along x=63.5mm, z=5mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=63.5, z=5, next line: zone t="data along x=63.5mm, z=10mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=63.5, z=10, next line: zone t="data along x=63.5mm, z=20mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
x=63.5, z=20, next line: zone t="data along z=0.4mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
z=0.4, y=0, next line: zone t="data along z=10mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data along 
z= 1 0, y=0 

b. Phase-averaged quantities at 9 different phases during the cycle: 0 deg, 40 deg, 80 deg, 120 
deg, 160 deg, 200 deg, 240 deg, 280 deg, 320 deg.; where you should align the phases of your 
computation as described below. Submit the following phase-averaged <> quantities: u/Uinf, 
v/Uinf, w/Uinf, u'u'bar/Uinf A 2, v'v'bar/Uinf A 2, w’w'bar/Uinf A 2, u'v'bar/Uinf A 2, u'w'bar/Uinf A 2, 
v'w’bar/Uinf A 2 (nondimensionalized), where: u = phase-averaged downstream velocity 
component, v = phase-averaged spanwise velocity component, w = phase-averaged vertical 
velocity component, u'u’bar = phase-averaged turbulent normal stress in downstream direction 
v'v'bar = phase- averaged turbulent normal stress in spanwise direction, w'w'bar = phase- averaged 
turbulent normal stress in vertical direction, u'v'bar = phase-averaged turbulent shear stress in x-y 
plane, u'w'bar = phase-averaged turbulent shear stress in x-z plane, v'w'bar = phase-averaged 
turbulent shear stress in y-z plane The locations for these data are the same as for the long-time- 
averaged quantities. Name these files: 

case2.phase000.ANYTHING.dat case2.phase040.ANYTHING.dat 

case2.phase080.ANYTHING.dat case2.phasel20.ANYTHING.dat 

case2.phasel60.ANYTHING.dat case2.phase200.ANYTHING.dat 

case2.phase240.ANYTHING.dat case2.phase280.ANYTHING.dat 

case2.phase320.ANYTHING.dat 

-where "ANYTHING" can be any descriptor you choose (should be different for each file if 
you are submitting multiple runs)-the file should be in 12-column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 
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3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #number of time steps per cycle (pound sign needed) 

6th line: #brief description of code/method (pound sign needed) 

7th line: #other info about the case, such as spatial accuracy (pound sign needed) 

8th line: #other info about the case, such as turb model (pound sign needed) 

9th line: #other info about the case (pound sign needed) 

10th line variables=:”x, mm”,”y, mm”, “z, mm”, "u/Uinf", "v/Uinf", "w/Uinf", 
"uu/Uinf", "w/Uinf', "ww/Uinf", "uv/Uinf", "uw/Uinf", "vw/Uinf", 

11th line: zone t="data along x=0, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 

ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=0, y=0 
next line: zone t="data along x=44.45mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 

ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=44.45, y=0 
next line: zone t="data along x=50.8mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=50.8, y=0 
next line: zone t="data along x=57.15mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=57. 15, y=0 
next line: zone t="data along x=63.5mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=63.5, y=0 
next line: zone t="data along x=76.2mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=76.2, y=0 
next line: zone t="data along x=101.6mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x= 10 1.6, y=0 
next line: zone t="data along x=57. 15mm, z=5mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=57. 15, z=5 
next line: zone t="data along x=57.15mm, z=10mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=57. 15, z=10 
next line: zone t="data along x=57.15mm, z=20mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=57. 15, z=20 
next line: zone t="data along x=63.5mm, z=5mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=63.5, z=5 
next line: zone t="data along x=63.5mm, z=10mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=63.5, z=10 
next line: zone t="data along x=63.5mm, z=20mm" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along x=63.5, z=20 
next line: zone t="data along z=0.4mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along z=0.4, y=0 
next line: zone t="data along z=10mm, y=0" 

subsequent lines: x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf uu/Uinf A 2 vv/Uinf A 2 
ww/Uinf A 2 uv/Uinf A 2 uw/Uinf A 2 vw/Uinf A 2 <- this is the data along z=10, y=0 
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The sample datafile case2.phase000.SAMPLE.dat can be downloaded from the website 

c. Phase-averaged time-history values of <u>, <v>, and <w> (nondimensionalized by Uinf) as 
a function of phase (deg) at three approximate point locations: (x,y,z)=(50.63, 0, 0.40)mm, 
(57.15, 0, 10)mm, and (63.5, 0, 10)mm. Give the data at every time step taken. In other words, 
if the time step yields 100 steps per cycle, then give 100 phases between 0 deg and 360 deg. The 
phases of your computation should be aligned as described below. Name this file: 
case2.phasehist.ANYTHING.dat-where "ANYTHING" can be any descriptor you choose 
(should be different for each file if you are submitting multiple runs)-the file should be in 7- 
column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #number of time steps per cycle (pound sign needed) 

6th line: #brief description of code/method (pound sign needed) 

7th line: #other info about the case, such as spatial accuracy (pound sign needed) 

8th line: #other info about the case, such as turb model (pound sign needed) 

9th line: #other info about the case (pound sign needed) 

10th line: variables= "phase, deg","x, mm","y, mm","z, mm","u/Uinf","v/Uinf","w/Uinf" 
11th line: zone t="x=50.63mm, y=0mm, z=0.4mm" 

subsequent lines: phase x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data at 
x=50.63, y=0, z=0.4 

next line: zone t="x=57.15mm, y=0mm, z=10mm" 

subsequent lines: phase x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data at 
x=57.15, y=0, z= 1 0 

next line: zone t="x=63.5mm, y=()mm, z=10mm" 

subsequent lines: phase x(mm) y(mm) z(mm) u/Uinf v/Uinf w/Uinf <- this is the data at 
x=63.5, y=0, z= 1 0 

The sample datafile case2.phasehist.SAMPLE.dat should be downloaded from the website 


d. Field line-contour-plots (in one of the following formats: ps, eps, or jpg) of long-time- 

averaged streamwise velocity (u/Uinf) in the planes y=0mm (centerplane), x=57.15mm (ID 
downstream), and x=76.2mm (4D downstream). These plots should be black-and-white line 
plots. The y=0 plane plot should go from approx x=45mm to 125mm and z=0 to 36mm. The 
x=const-plane plots should go from approx y=-18mm to +18mm, and z=0 to 36mm. The x-to-z 
ratio of the plots should be 1.0. Plot u/Uinf line contours of -0.5 through 1.5 in increments of 0.1. 
Label the contour lines, if possible. The purpose of submitting these plots is to get a qualitative 
picture of the long-time-averaged flowfield, indicative of the dynamic range of the orifice's 
influence. Altogether, submit 3 plot files. Name these files: case2.uavg.yO.ANYTHING.eps 
case2.uavg.x57.15.ANYTHING.eps case2.uavg.x76.2. ANYTHING.eps (where the 

"eps" in this case means encapsulated postscript - use ps, or jpg instead if appropriate). 

e. Field line-contour-plots (in one of the following formats: ps, eps, or jpg) of phase-averaged 

streamwise velocity (<u>/Uinf) in the planes y=0mm (centerplane), x=57.15mm (ID 
downstream), and x=76.2mm (4D downstream) at the following phases: 40 deg, 120 deg, 200 
deg, and 280 deg.; where you should align the phases of your computation as described below. 
These plots should be black-and-white line plots. The y=0 plane plots should go from approx 
x=45mm to 125mm and z=0 to 36mm. The x=const-plane plots should go from approx y=- 
18mm to +18mm, and z=0 to 36mm. The x-to-z or x-to-y ratio of the plot should be 1.0. For all 
plot files, plot <u>/Uinf line contours of -1.0 through 2.0 in increments of 0.1. Either (a) plot the 
lines of negative velocity as dashed lines and the lines of positive velocity as solid lines, or (b) 
label the contour lines, or (c) do both. The purpose of submitting these plots is to get a 
qualitative picture of the phase-averaged flowfield at particular selected times of interest. 
Altogether, submit 12 plots files. Name these files: case2.phase040.y0.ANYTHING.eps 
case2.phasel20.y0.ANYTHING.eps case2.phase200.y0.ANYTHING.eps 
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case2.phase280.y0.ANYTHING.eps case2.phase040.x57.15.ANYTHING.eps 

case2.phasel20.x57.15.ANYTHING.eps case2.phase200.x57.15.ANYTHING.eps 

case2.phase280.x57.15.ANYTHING.eps case2.phase040.x76.2. ANYTHING.eps 

case2.phasel20.x76.2. ANYTHING.eps case2.phase200.x76.2. ANYTHING.eps 

case2.phase280.x76.2.ANYTHING.eps (where the "eps" in this case means encapsulated 
postscript - use ps, or jpg instead if appropriate). 

Definition of Phase for the Computations 

Matching the same phases with the experiment is not necessarily straightforward. One way to 
do it is to try to align a quantity from experiment (such as vertical velocity near the jet exit), but 
this can be imprecise because the CFD and experimental data are not necessarily well-behaved 
sine-waves. The best criteria for determining phase may in fact be a different measure altogether. 

On the other hand, in this workshop we want to be able to compare CFD results with each 
other, so it is important to try to achieve the same phase definitions in order that all computations 
are similarly aligned. We have tried to choose criteria for determining phase that approximates 
experiment AND is specific enough so that different CFD solutions can be meaningfully 
compared. 

Therefore, although participants are given some latitude to determine phase as appropriate, we 
encourage everyone to use the following steps to define phase in a uniform fashion for Case 2: 

• Step 1. Output vertical phase- averaged velocity (w) at the following point in space (over 
the orifice) as a function of your iteration or time step number: (x,y,z)=(50.63, 0, 0.40)mm. Find 
wmax and wmin over the course of one phase-averaged period. 

• Step 2. Compute the mid-value wavg=(wmax+wmin)/2 

• Step 3. Define Phase=50 deg as the time when your velocity at this location 
approximately equals wavg (INCREASING). All other phases can be referenced from this, via 
the following relation: 

Phase=(iter-it50)*360/nstep+50 

where: 

iter =iteration (or time step) number 

nstep=no. of time steps per cycle 

it50=iteration number when Phase=50 according to the above criteria 

For example, if you are running 360 steps per cycle and you match the above criteria at time 
step number 5492, then 

Phase=iter-5442 

Thus, when iter=5442, then Phase=0; when iter=5622, then Phase=180; when iter=5802, then 
Phase=360. Note that Phase=360 also corresponds with Phase=0 (it repeats every 360 deg). This 
is illustrated in the following figure (updated on 24 December 2003): 
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Example of defining phase from phase -averaged 



iteration number 


As another example, say the computation ran 1080 steps per cycle and you match the above 
criteria at time step number 10002, then 
Phase=(iter- 1 0002)/3+50 

Thus, when iter=9852, then Phase=0; when iter=10392, then Phase=180; when iter=10932, 
then Phase=360. 
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CASE 2: SYNTHETIC JET IN A CROSSFLOW 


Julien Dandois t , Eric Gamier * and Pierre Sagaut § 
ONERA, BP 72, 29 av. de la division Leclerc, 92322 Chatillon Cedex, France 


1 Introduction 

Three-dimensional, compressible simulations of the interaction between a synthetic jet and a turbulent bound- 
ary layer have been performed using Large Eddy Simulation. A modeled boundary condition has been used 
for the actuator. A study of the effects of the time step, number of sub-iterations in the time integration pro- 
cess, actuator boundary condition has been investigated. 

2 Solution Methodology 

The unsteady, three-dimensional, compressible Navier-Stokes equations are solved using Large-Eddy Sim- 
ulations. Any flow variable can be written as + <&', where represents the resolved part of the 

variable and <J>' its subgrid part. The filtering operator is classically defined as a convolution product on the 
grid cell. The FLU3M solver, developed by ONERA, is based on a cell centered finite volume technique 
and structured multiblock meshes. For efficiency, an implicit time integration is employed to deal with the 
very small grid size encountered near the wall. The time integration is carried out by means of the backward 
scheme of Gear which is second-order- accurate: 

dw — 2 w n + \w n ~ l 

dt st y J 

Because of the implicit terms, a non-linear system has to be solved. The Newton-Raphson method is used 
to compute w nJrl . At each iteration of this inner process, the inversion of the linear system relies on lower- 
upper symmetric Gauss-Seidel implicit method. More details about these numerical points are available in 
Ref. [3]. A special effort has been carried out to minimize the intrinsic dissipation of the scheme and the 
computational cost. The spatial scheme is the scheme proposed by Mary and Sagaut [2]. It is based on the 
AUSM+(P) [4] scheme, whose dissipation is proportional to the local fluid velocity so it is well adapted to 
low-Mach number boundary-layer simulations. As we are not interested by the shock capturing properties 
of the scheme, simplified formulas have been developed to approximate the Euler fluxes: 

F[ +h = U^Ql + Q R )/2 - \U dts \(Q R - Q L )/2 + P (2) 

where U% denotes the interface fluid velocity, L[ R the left and right third-order MUSCL interpolation. The 
state vector Q is defined as (p, pu \ , pu 2 , pu 3 , pE + p) 1 , whereas the pressure term P is given by [0, (pl + 
Pr)/ 2, 0,0,0]. The symbol u^s, which indicates a local fluid velocity, characterizes the numerical dissipa- 
tion acting on the velocity components. It is defined as follows: 

Udis = rnax(\ui L + u3r\/2, c x ) (3) 

tphD student, ONERA, Applied Aerodynamic Department 
+ Research Engineer, ONERA, Applied Aerodynamic Department 

§ Professor, Laboratoire de Modelisation pour la Mecanique and ONERA, CFD and Aeroacoustics Department 


2.3.1 



where c\ is a constant parameter. To enforce the pressure/velocity coupling in low-Mach-number zone, a 
pressure stabilization term is added to the interface fluid velocity: 

U\ = {uil + uir)/2 - c 2 (pr - pl) (4) 

where c 2 is a constant parameter. For an accuracy reason the values of ci and c 2 should be chosen as small 
as possible to minimize the numerical dissipation. For a stability reason these parameters cannot be smaller 
than 0.04. For the viscous fluxes, a second-order-accurate centered scheme is used. 

3 Model Description 

The subgrid scale model used is the selective mixed scales model [1]. This model has been shown to be 
effective in turbulent flows. The expression of the subgrid viscosity is: 

Vsm = C m (a)\S(~X,t)\ a 2 A 1+a f s (5) 

where C m = 0.0, o = 0.5, f s is a selection function which tests the tridimensionality of the flow (see [1] for 
more details) and q c is the kinetic energy of the smallest resolved scales which is evaluated with a test filter 
with a cutoff A = \/6 A where A is the cubic root of the cell volume: 

^ (u — u) (u — u) where Ui = ^~ 1 ^ ^t+i (6) 

4 Implementation and Case Specific Details 

Table 1 provides the mesh size and grid spacings in the central refined zone of the inflow domain. 


Table 1 : Computational mesh parameters 


Ax (mm) 

Ay (mm) 

^%min (mm) 

Ax + 

Ay+ 

A^ + 

1.16 

0.44 

0.036 

100 

36 

3 


The size of the computational domain is the same than the one of the RANS grid posted on the workshop 
website: the streamwise length is 152.4mm or 24D (8D upstream and 16D downstream), the spanwise and 
vertical length are 76mm or approximately 12D. There are 180 points on the jet circumference. The grid is 
composed of 14 domains. The cells number repartition is shown in Table 2. 


Table 2: Cells repartition 


External field 

Orifice Cavity Total 

1.25 million 

263 000 98 000 1.7 million 


Figures 1 and 2 provide a view of the grid in the x-y and x-z planes. 





Figure 1: X-Y view of the grid. 



Figure 2: X-Z view of the grid. 


The boundary conditions are shown in Fig. 3. On the whole cavity’s bottom surface, a bio wing/ suction 
condition with a top-hat distribution which varies sinusoidaly in time is implemented in order to simulate the 
diaphragm movement: u(x,t) — ( o * cos(27rft) with Uq = 0.273m.s -1 and / = 150 Hz. f/ 0 has been 
calculated using, on the one hand, the sections ratio between the orifice and the cavity bottom surface and, 
on the other hand, a target output velocity of 50 m.s ~ l . 


subsonic outflow 



blowing/suction 


Figure 3: Boundary conditions. 




A study of the influence of the actuator boundary condition, time step and number of sub-iterations has 
been performed. For those case, the blowing/suction condition is applied both on piston and diaphragm. For 
the case D, instead of using a sinusoidal function, a polynomial fit of driver displacement derivative given 
on the workshop’s website has been used. In that case, the blowing/suction condition is applied only on the 
piston. The inflow boundary condition is based on a steady RANS mean velocity profile. Reaslistic turbulent 
inflow conditions will be implemented in a follow-up of this study. The table 3 shows the different varying 
parameters. 


Table 3: Computational parameters of the study 


Case designation 

Time step 

Number of 
subiterations 

Surface of the 
boundary condition 

Uq ( m.s 

A 

1 flS 

4 

piston+diaphragm 

0.273 

B 

1 jlS 

6 

piston+diaphragm 

0.273 

C 

0.5 jis 

4 

piston+diaphragm 

0.273 

D 

1 flS 

4 

piston 

0.725 


The figure 4 shows that the solution is not very sensitive to the investigated parameters. For the submit- 
ted data, the time step was fixed to 1 jis and the number of sub-iterations has been increased to 8 because 
residues, not shown here, were not sufficiently converged with 6 sub-iterations. On the plot of figures 4 and 
5, the phase duration is 10 deg like in the experimental results provided in the website. 



Figure 4: Influence of the time step and the number of sub-iterations on the phase averaged values of W/Ui n j 
after 2 periods at (y=0, z=0.4mm). 


Although the actuator boundary condition of the case D is more realistic, the results are less accurate than 
those of the case A (see fig. 5) and the simpler boundary conditions has been finally selected. 



Figure 5: Influence of the actuator boundary condition on the phase averaged values of IF/ ( ' n/ j after 2 pe- 
riods at (y=0, z=0.4mm). 


For the definition of phase, the output vertical velocity has been plotted as a function of time (see fig. 6). 
The definition given on the workshop website has been used. The submitted data have been phase-averaged 
on 8 periods. 



Figure 6: Vertical velocity just above the orifice upstream of the center (x,y,z)=(50.63,0,0.4)mm. 
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CASE 2: URANS APPLICATION WITH CFL3D 


C. L. Rumsey 


Computational Modeling & Simulation Branch, NASA Langley Research Center, Hampton, VA 23681-2199 


Solution Methodology 

This case was run using CFL3D, a multi-zone Reynolds-averaged Navier-Stokes code developed at NASA 
Langley [1]. It solves the thin-layer form of the Navier-Stokes equations in each of the (selected) coordinate 
directions. It can use 1-to-l, patched, or overset grids, and employs local time step scaling, grid sequencing, 
and multigrid to accelerate convergence to steady state. In time-accurate mode, CFL3D has the option to 
employ dual-time stepping with subiterations and multigrid, and it achieves second order temporal accuracy. 

CFL3D is a finite volume method. It uses third-order upwind-biased spatial differencing on the con- 
vective and pressure terms, and second-order differencing on the viscous terms; it is globally second-order 
spatially accurate. The flux difference- splitting (FDS) method of Roe is employed to obtain fluxes at the 
cell faces. It is advanced in time with an implicit three-factor approximate factorization method. 


Model Description 

For this test case, two different turbulence models have been run to date. The first is the one-equation 
Spalart-Allmaras model (SA) [2], and the second is the two-equation shear- stress transport model of Menter 
(SST) [3, 4]. These are both linear eddy-viscosity models that make use of the Boussinesq eddy-viscosity 
hypothesis. The equations describing these two models can be found in their respective references. 

In CFL3D, the models are implemented uncoupled from the mean-flow equations. They are solved using 
a three-factor implicit approximate factorization approach. The advection terms are discretized with first- 
order upwind differencing. The production source term is solved explicitly, while the advection, destruction, 
and diffusion terms are treated implicitly. 


Implementation and Case Specific Details 

In this flow, a synthetic jet issues into a turbulent boundary layer through a circular orifice of diameter 6.35 
mm in the floor. The flow is characterized by a forcing frequency of 150 Hz, with a maximum discharge 
vertical velocity of approximately 1 .25/ x . The approach boundary layer thickness is somewhat greater than 
20 mm, the freestream Mach number is M — 0.1, and the Reynolds number is approximately Re = 14160 
per orifice diameter. 

The grid used was the supplied structured grid number 1 (which contains 7-zones connected in a 1-to-l 
fashion, and approximately 4.09 million grid points), as well as a medium-level grid made from the fine grid 
by extracting every-other point in each coordinate direction. The SA model was solved on both the fine and 
medium grids, whereas the SST model was only solved on the medium grid. 

The time step chosen was one that yielded 720 time steps per cycle of the forcing frequency. For the SA 
model, 5 subiterations were employed per time step. For the SST model, however, a large instability was 
noted when 5 subiterations were employed. This instability showed up as a series of un-physical oscillations 
in the flowfield variables at given points in space during part of the cycle. By increasing the number of 
subiterations to 10, the level of this instability decreased, but it did not go away. See for example, Fig. 1, 
which shows the vertical (W) velocity component as a function of phase ID downstream of the center of the 
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Figure 1: Effect of subiterations on result using SST at x=57.15 mm, z=10 mm (ID downstream), medium 
grid. 


orifice, at z=10 mm above the floor. Therefore, the current SST results should be considered as incomplete; 
future work will address eliminating the instability. 

The boundary conditions were as follows. At the tunnel floor as well as on the inner lip of the orifice 
and upper wall of the cavity, solid wall adiabatic boundary conditions were applied. The side walls of 
the cavity were modeled as slip walls, and the bottom surfaces of the cavity employed a time-dependent 
boundary condition. In other words, the moving diaphragm was modeled on the stationary grid through the 
use of a boundary condition that imposed a vertical velocity in a sinusoidal manner. (Note that the bottom 
of the cavity in the supplied grid modeled both the elastic diaphragm as well as the solid plate glued to its 
center. The time-dependent boundary condition was applied both on the part modeling the elastic diaphragm 
as well as on the part modeling the solid plate.) The time-dependent boundary condition set the velocity 
components as follows: 


[7 = 0 V = 0 W = [(pW) ma Jp\cos(27rFt) (1) 

where F is the frequency and t is the time. With this boundary condition, the density and pressure are 
extrapolated from the interior of the domain. The (pW) max was chosen by trial and error in order to achieve 
an approximate match of the vertical velocity component at the outflow of the orifice with the experiment 
(the final value used was (pW) max = 0. 0008 p^a^ where a ^ is the reference speed of sound). The 
resulting W and U components of velocity at the orifice exit can be seen in Figs. 2 and 3. Note that use 
of the current boundary condition fails to capture the sharp peak in W and the additional hump on the 
downstroke. Also, it should be noted that the experimental data exhibited significant side (V -component) 
velocities during the cycle at this location, whose cause has not been accounted for. The current CFD method 
makes no attempt to duplicate this V -velocity component. Also shown in these figures is the effect of fine 
vs. medium grids, which is very small at this point in the flowfield. 

The supplied grid extends .0508 m upstream of center of the orifice, and .1016 m downstream of the 
center of the orifice. Its height is .076 m above the floor, and its width extends from y=-.038 m to +.038 m 
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Figure 4: Contours of U/U^ at x=76.2 mm (4D downstream), phase=200°, SA on fine grid. 


(total width of .076 m). Note that this height and width are smaller than the wind tunnel height (approx .249 
m) and width (.381 m). The boundary conditions at the top, side, and downstream faces of the grid in the 
tunnel are a farfield Riemann-type. At the upstream face in the tunnel, the density and velocity components 
are specified in order to approximately match the experimental boundary layer thickness, a nd th e turbulence 
data is specified in order to approximate a fully-developed turbulent boundary layer in its u f w f component. 
The pressure at the inflow is extrapolated from the interior of the domain. 

Examples of the effects of grid and turbulence model can be seen in Figs. 4, 5, and 6, at phase=200° in 
the plane 4D downstream (x=76.2 mm). SA on the fine grid yields a somewhat rounder structure than SA on 
the medium grid, but overall the results are very similar. The differences between the SA and SST models 
(on the medium grid) are also relatively minor, with the SST model producing a thinner and taller structure. 
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CASE 2: SYNTHETIC JET IN A CROSSFLOW 


A. Azzi 1 , and D. Lakehal 2 


l LaboratoiredeMecaniqueAppliquee, FacultedeGenie-Mecanique, Universite des Sciences et de la 
Technologie d’Oran, USTO, BP 1505 El-Mnaouar, Oran, Algeria 
2 Institut of Energy Technology, ETH Zentrum/CLT2, CH-8092 Zurich, Switzerland. 


Solution Methodology 

The governing Unsteady Reynolds Averaged Navier Stokes (U-RANS) equations are solved by the use of 
a three-dimensional finite-volume method that allows the use of arbitrary nonorthogonal grids, employing 
a cell-centered grid arrangement. A detailed description of the method is reported in Majumdar et al. [1], 
and the multi -block technique which was introduced afterwards, in Lakehal et al. [2]. The momentum- 
interpolation technique is used to prevent pressure-field oscillations tending to appear in the cell-centered 
grid arrangement. The pressure-velocity coupling is achieved by using the well-known SIMPLEC 
algorithm. The present computations were performed employing the second-order Centered Differential 
Scheme for all variables applied by a deferred-correction procedure. The resulting system of difference 
equations was solved using the strongly implicit procedure (SIP) algorithm. The unsteady terms are 
discretized by an implicit second order time scheme and computations are conducted until a nearly 
repeating periodic solution is reached. Finally, the long-time-averaged (time-independent mean) values 
for each quantity are computed as an average over the last period of oscillation. 


Model Description 

The Reynolds-stress tensor is approximated within the context of the k — £ turbulence model, 
considering both its linear and nonlinear forms. For the sake of consistency, the two models are coupled 
with a one-equation model resolving the near-wall viscosity affected regions. 

The Two-Layer DNS-Based k — £ Turbulence Model (TLV) 

The two-layer approach represents an intermediate modeling strategy between wall function and pure 
low-Re number models. It consists in resolving the viscosity-affected regions close to walls with a one- 
equation model, while the outer core flow is treated with the standard k — £ model. 

The two-layer model used here is a re-formulation of the so-called velocity-scale-based model (TLV) of 
Rodi et al. [3], in the sense that the turbulent kinetic energy is re-incorporated as a velocity scale. A 
detailed description of the model is reported in Azzi and Lakehal [4]. 

The EASM model (GSLT) 

In this model, the Reynolds stress is represented algebraically in terms of a series of combinations 
between vorticity and strain. With this, the model is expected to reproduce secondary flows with better 
accuracy. Although these relationships may be obtained by invoking various strategies, their common 
starting point is the assumption of homogenous turbulent flows in the limit of equilibrium. 

The model used in the present study is a modified version of the EASM model of Gatski and Speziale [5]. 
In a recent effort to make the model applicable to broad range of practical flows, Lakehal and Thiele [6] 
emphasized two important aspects: (i) the development of a generalized relation for P k j£ , and (ii)the 
formulation of a better regularization procedure for C u . 
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Implementation and Case Specific Details 

The computational domain is composed by three blocks, including the jet hole and the internal cavity. 
Only half of the physical domain is computed and a symmetry boundary is applied at y=0 plane. The 
computational domain extends from -8 diameters upstream to +20 diameters downstream of the centre of 
the hole. It extends up to 8 diameters over the flat plate in the vertical direction and up to 4 diameters in 
the span wise direction (figure 1). A test with 6 diameter in span wise direction was also tested but gives 
similar values as with 4d. 

Conforming to the guideline instructions, the coordinate system is set with x downstream, z up, and y 
span wise, with the (x,y,z)=(0,0,0) origin at 8 diameters (50.8 mm) directly upstream of the center of the 
orifice (the orifice diameter is = 6.35 mm). 

Results are provided as long-time averaged and phase averaged quantities along specified lines as showed 
in figure 2. 

Boundary conditions: 

The moving diaphragm is represented by an unsteady boundary condition imposed at the bottom side of 
the cavity. The averaged injected velocity is computed using a blowing ratio of nearly 1, and the unsteady 
condition is imposed using 

V=Vcos(2jt.150 .t) (1) 


where 150 is the diaphragm frequency. 


Preliminary computations have been performed and velocities components at point (50.63 , 0.0 , 0.4) have 
been found to compare well the averaged data (Figure 3). 

The inlet boundary conditions provided in the file (UpstreamBC.dat) are used as mainstream flow inlet 
conditions. The kinetic turbulent energy is computed according the Reynolds stress components given in 

/ du^ 

(UpstreamBC.dat) file and its rate of dissipation is computed assuming P = £ and £ = u V . 

A ratio of turbulent to molecular viscosity of 20 is used for injected jet assuming a turbulence intensity of 
2%. The time step is fixed at dt =1.852E-5, which gives 360 steps per cycle. 
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Figure 2: Lines where results are provided. 






Figure 3: Streamwise and vertical velocities just above the orifice, 
0.17 upstream of its centre : (x,y,z)=(50.63, 0, 0.4) mm 
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CASE 2: SIMULATION OF A PERIODIC JET IN A CROSS FLOW 
WITH A RANS SOLVER USING AN UNSTRUCTURED GRID 


H. L. Atkins 
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Introduction 

A second-order unstructured-grid code, developed and used primarily for steady aerodynamic simulations, is 
applied to the synthetic jet in a cross flow. The code, FUN3D, is a vertex-centered finite-volume method orig- 
inally developed by Anderson[l, 2], and is currently supported by members of the Fast Adaptive Aerospace 
Tools team at NASA Langley. Used primarily for design[3] and analysis[4] of steady aerodynamic configu- 
rations, FUN3D incorporates a discrete adjoint capability, and supports parallel computations using MPI. 


Solution Methodology 

A detailed description of the FUN3D code can be found in the references given above. The code is under 
continuous development and contains a variety of flux splitting algorithms for the inviscid terms, two meth- 
ods for computing gradients, several turbulence models, and several solution methodologies; all in varying 
states of development. Only the most robust and reliable components, based on experiences with steady 
aerodynamic simulations, were employed in this work. 

As applied in this work, FUN3D solves the Reynolds averaged Navier-Stokes equations using the one 
equation turbulence model of Spalart and Allmaras[5]. The spatial discretization is formed on unstructured 
meshes using a vertex-centered approach. The inviscid terms are evaluated by a flux-difference splitting 
formulation using least-squares reconstruction and Roe-type approximate Riemann fluxes. Green-Gauss 
gradient evaluations are used for viscous and turbulence modeling terms. 

The discrete spatial operator is combined with a backward time operator which is then solved iteratively 
using point or line Gauss-Seidel and local time stepping in a pseudo time. For steady flows, the physical 
time step is set to infinity and the pseudo time step is ramped up with the iteration count. A second-order 
backward in time operator is used for time accurate flows with 20 to 50 steps in the pseudo time applied at 
each physical time step. 

For this effort, FUN3D was modified to support spatially varying boundary and initial conditions, and 
unsteady boundary conditions. Also, a specialized in/out flow boundary condition was implemented to 
model the action of the diaphragm. This boundary condition is described below in more detail. 

The grids were generated using the internally developed codes GridEX[6] for meshing the surfaces and 
inviscid regions of the domain, and for CAD access; and MesherX[7] for meshing the viscous regions. Grid 
spacing in on the surfaces and in the inviscid regions are indirectly controlled by specifying sources. The 
viscous layers are generated using an advancing layer technique. MeshersX allows the user to control the 
spatial variation of the first step off the surface, growth rates, and the termination criterion by providing 
small problem dependent subroutines. 


Modeling of Diaphragm Boundary 

A specialized in/out flow boundary condition is formulated to model the action of the diaphragm at a station- 
ary boundary. This boundary condition is similar to a standard characteristic boundary condition, commonly 


2 . 6.1 



applied at far- field boundaries, in that the inviscid boundary flux is evaluated from an intermediate solution 
state Ui that is computed from the current boundary solution Ub and a prescribed external state U v . In the 
weak form implemented in FUN3D, the boundary solution is not replaced by the intermediate solution but is 
allowed to evolve. The viscous gradients are evaluated using only the boundary and interior solution values, 
and have no knowledge of the prescribed or intermediate solution states. 

At a far-held boundary, the intermediate solution is determined by the characteristic relationships that 
model the convective and acoustic waves that are assumed to exist between the boundary and the far-held. 
For example, at a subsonic boundary 


= n-m = irm and iw.> = {2$|i ^“° «> 

where R + (~) denotes the characteristic variable associated with acoustic waves leaving(entering) the do- 
main, R c denotes characteristic variables associated with convective waves, and V n is the normal velocity 
through the boundary, with outhow taken as positive. 

In the case of a moving diaphragm, the “prescribed” state is that at the diaphragm face, and only the 
diaphragm velocity is known. If the boundary were moving with the diaphragm, than the intermediate state 
would be given by 


v n = 0, R + (Ui) = R + {U b ), and R c {U t ) = R c {U b ) . (2) 

Simply applying these conditions at a stationary location, by setting V n equal to the diaphragm velocity, 
would result in under specifying the flow during the inflow phase of the simulation. To stabilize the simu- 
lation, the latest value of R c is saved during the outflow phase of the calculation, and reapplied during the 
inflow phase to give: 

Vn = acosfojt), = and R C (U,) = { ^=° (3) 

where R c denotes the value saved from the most recent outflow cycle. In practice, only the entropy is saved 
and reapplied. The boundary tangential velocity is allowed to evolve without constraint. Also, during the 
inflow cycle, the entropy is gradually relaxed back toward its initial value. 


Implementation and Case Specific Details 

The geometry of the cavity is simplified by making the diaphragm flat with the height of the diaphragm 
chosen so that the “at rest” volume of the cavity is unchanged. To reduce the problem size, the domain is 
divided at the tunnel centerline and only half of the domain is grided. Although all simulations presented 
here were performed on this half-domain geometry, it is possible to reflect the grid about the centerline to 
obtain a symmetric grid for the larger problem. The computational domain extends from 8 jet diameters 
upstream and to the side of the jet center, and to 16 jet diameters downstream of the jet center. 

The fine grid was sized so that the first spacing on the tunnel wall would have a y+ < 1. Spacings 
inside the cavity were based on preliminary simulations with steady blowing applied at the diaphragm face. 
Griding sources were placed around the lip of the jet and along the jet centerline in an effort to improve 
the clustering there. Figures l(a-c) shows three views of the mesh on the symmetry plane that illustrate the 
mesh distribution near the jet. 

A coarse grid was generated simply by doubling all grid sources and modifying the growth rate of 
the viscous layers. Originally, this grid was intended only to provide rapid turn-around while sorting out 
boundary and initial conditions. The fine grid has 1457853 tetrahedra and 255426 nodes. The coarse grid 
has 254046 tetrahedra and 46063 nodes. 

Characteristic boundary conditions are applied weakly at the inflow, outflow and top boundaries. Sym- 
metry conditions are imposed at both y=constant boundaries. No-slip conditions are enforced on the tunnel 


2 . 6.2 





Figure 1: Grid on the symmetry plane: a) overall view, b) near jet exit, c) near corner of jet exit. 


floor, the top wall of the cavity, and in the jet contraction. At these no-slip boundaries, the temperature is 
set to the adiabatic wall temperature, and the density is determined from the continuity equation. The side 
walls of the cavity are treated as inviscid walls, which FUN3D enforces weakly. The in/out flow boundary 
condition described in the previous section is applied on the lower wall of the cavity. 

The inflow and initial conditions in the tunnel region were generated by performing a separate boundary 
layer simulation with FUN3D and extracting the solution at the appropriate Reynolds number. Initially, there 
is no flow in the cavity and the pressure is set to the freestream value. The simulation was started impulsively 
with the diaphragm velocity (normalized by u^) specified as -0.007 cos(u;£). The simulations used 720 
time steps per period, with results saved every 5 degrees of phase. The long time averages requested by the 
workshop were computed by averaging these 5-degree samples. The fine grid simulations required 35 hours 
per period when using 16 intel processors; the coarse grid simulations required less than 1 1 hours per period 
when using 8 processors. 

The fine and coarse grid simulations produced noticeably different solutions in several respects. Fig- 
ures 2 (a) and (b) show the fine grid solution history at a location 0.17mm upstream of the jet center. 
Although the vertical velocity component W settles quickly to a periodic solution, the streamwise velocity 
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Figure 2: Solution at 0.17mm upstream of jet center: a) and b) fine grid startup history of W and U, 
respectively; c) comparison of U. 


component U requires 4 periods. In the coarse grid simulation, the solution requires even longer, 6 periods, 
before the U velocity component becames periodic. The fine and coarse grid simulations give similar results 
for W, but the U velocity components, compared in fig. 2(c), have completely different character during 
the blowing phase. The coarse grid produces results similar to the experiment, while the fine grid predicts 
a large negative stream wise velocity component during the blowing cycle. The contour plots of U, shown 
in fig. 3, indicate that this negative streamwise velocity is due to a vortical behavior that develops in the jet 
exit flow. 
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Figure 3: Streamwise velocity contours from fine grid simulation during peek blowing: phase = 160. 
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Introduction 

Active flow control technologies is growing area of aerodynamic research in the early 21 st 
century. The goal is to prevent boundary layer separation and as such it is often applied to designs 
of high-lift airfoils. Forced oscillations superimposed on a mean flow which is on the verge of 
separation point is an effective means to delay boundary-layer separation, such as blowing or 
suction techniques. However, the progress in active flow control technologies has often been 
paced by the development of actuator capabilities. A popular current actuator is the synthetic jet, 
which has demonstrated capabilities regarding separation control and thrust vectoring. 

Over the past several decades, both computational fluid dynamics (CFD) algorithms and 
computer technologies have progressed tremendously. These advances have made unstructured 
grids more attractive with their ability to smoothly conform to varying flow conditions and 
complicated boundaries with a single grid. However, challenges remain, including grid generation 
for a computational domain with complex geometries, well-balanced grid decomposition on 
distributed system, and efficient parallel performance. In this paper, a 3D numerical simulation of 
a synthetic jet into a cross flow using a new CFD code called UNCLE is described. UNCLE is a 
2D/3D finite volume unstructured unsteady incompressible Navier-Stokes solver. In order to take 
care of turbulence flow in most realistic cases, F. R. Menter’s shear-stress transport (SST) 
turbulence model [1] is employed to UNCLE. It is designed to study the challenges of using 
unstructured grid codes on high-performance parallel computers to simulate realistic cases. To 
increase flexibility in complex geometries, UNCLE may use a variety of grid types, such as 
triangles, quadrilateral, tetrahedron and hexahedra meshes. In order to achieve good load balance 
for parallel computing, METIS [2] is used to partition the grid. 

Solution Methodology 


UNCLE employs a pressure-based SIMPLE algorithm with second order accuracy in both time 
and space. A second order upwind scheme is used for computing advection terms. Non-staggered 
grids with the Rhie and Chow momentum interpolation method [3] are employed to avoid 
checkerboard solutions. In order to take care of turbulence flow in most realistic cases, F. R. 
Menter’s shear-stress transport (SST) turbulence model [1] is implemented. It is designed to study 
the challenges of using unstructured grid codes on high-performance parallel computers to 
simulate realistic cases. To increase flexibility in complex geometries, UNCLE may use a variety 
of grid types, such as triangles, quadrilateral, tetrahedron and hexahedra meshes. 

In order to achieve good load balance for parallel computing, METIS [2] is used to 
partition the grid. Generally, there are two different partitioning approaches - vertex based and 
element based partitioning for mesh-partitioning as shown in Fig. 1(a) and (b) respectively. For 
vertex-based partitioning, the boundary elements are doubled and the vertices at the boundary are 
overlapped. Since the control volumes at boundary are not partitioned, only communication of 
boundary nodal properties is required. For element-based partitioning, the boundary vertices are 
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doubled. Because the control volumes at the boundary are split, all nodal points surround a 
boundary vertex are needed to interpolate the properties of the boundary vertices. Communication 
of boundary nodal properties in element based partitioning is heavier than vertex-based 
partitioning. On the other hand, vertex-based partitioning has to handle doubled elements at the 
boundary, it still costs computational time. For the purpose of direct use the information from 
METIS, element based partitioning is used for pre-processing code of UNCLE. 

Model Description 

F. R. Menter’s SST turbulence model [1] is employed in this simulation. 


Implementation and Case Specific Details 

In this case, the orifice diameter, 6.35 mm is chosen to be the reference length. Air density at sea 
level is 1.185kg/m 3 , viscosity is 18.4e-6kg/m-s, and the reference velocity is the freestream 
velocity which is 34.6m/s. According to the information above, the Reynolds number for this case 
is 1414.8. 

We assume the flow is symmetric about y = 0 plane, so only half of the experimental 
domain is used as the computational domain. An unstructured grid is generated to fulfill the grid 
format of UNCLE by using GAMBIT. The total grid has 0.3 million cells. The geometry of the 
grid is similar to the provided grids from the cfdval2004 website 

In order to simulate the moving diaphragm, a periodic velocity boundary condition is 
used to replace moving boundary condition. The non-dimensional period T can be calculated by: 

l/T = f* =f*£/u (1) 

where / is non-dimensional frequency, frequency / =150Hz, reference length /? =0.00635 m, 
and reference velocity u =34.6 m/s. As a result of eq. (1), the non-dimensional period T is 36.325 
approximately. Using the driver’s boundary conditions from the experimental data provided on 
the cfdval2004 website, the velocity at the center of the diaphragm is obtained. The time- 
dependent velocity at the center of the diaphragm can be approximated by curve-fitting method, 
and its mathematical formulation is described as: 

w(t) = a + b* cos (c *t + d), (2) 

a=2. 64228 le-6, b=0.020915027, c=0.17405397, d=-0.027378007 
where a, b, c, and d are dimensionless parameters and t is dimensionless time. Symmetry 
boundary condition is used along the symmetry plane. For the rest of the boundary, the no-slip 
condition is imposed as the wall boundary condition. In this simulation, the non-dimensional time 
for one step is chosen as 0.009 so that the total number of time steps for each cycle is 4036. 
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Fig. 1(a) Schematic of vertex based partitioning, (b) Schematic of element 
based partitioning. 
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Introduction 

The approach described here is to model the response of a synthetic jet to a grazing boundary layer using 
lumped element modeling. This effort complements more rigorous and expensive numerical simulations. 

Solution Methodology 

The approach used is to model the actuator orifice impedance with a grazing boundary layer and is based 
on a lumped element modeling (LEM) technique, following the recent paper by Gallas et al. [1]. 

Model Description 

For a full discussion on the LEM technique used herein, the reader is referred to the write-up for Case 1 
and to the details and references contained in Gallas et al. [1]. In summary, in LEM the individual 
components of a synthetic jet are modeled as elements of an equivalent electrical circuit using conjugate 
power variables (i.e., power = generalized “flow” x generalized “effort” variables); the lumped 
parameters represent generalized energy storage elements (i.e., capacitors and inductors) and dissipative 
elements (i.e., resistors), as shown in Figure 1. The frequency response function of the circuit is then 
derived to obtain an expression for Q out , the volume flow rate. LEM provides a compact nonlinear 
analytical model and valuable physical insight into the dependence of the device behavior on geometry 
and material properties. 



is the cavity acoustic impedance, Z a0 +Z aBL = AP out /Q out is the orifice acoustic impedance. 

Case 2 is similar to the device studied and modeled by Gallas et al., although now the driver is a voice- 
coil driven piston, and the synthesized jet interacts with a crossflow boundary layer. No information is 
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available about the dynamics of the electromagnetic transducer. Instead, the measured sinusoidal 
displacement of the (assumed) rigid piston is used to calculate the piston volume velocity 

Q d = v = j®Sd W 0 sm(cot), ( 1 ) 

where AV represents the volume displaced by the piston, S d is the effective moving area of the piston 
and W 0 is the amplitude of oscillation of the piston. This representation enables us to bypass the need of 
an expression for the impedance of the transducer shown in Figure 1 . 


The model considered and developed in [1] is for a device exhausting into a quiescent medium. Thus, in 
order to account for the incoming crossflow effect, a boundary layer impedance is introduced in series 
with the orifice impedance since they share the same volume flow Q out . The work done in the acoustic 
liners community has been used to derive a simple analytical expression for this grazing impedance. 
Specifically, it is derived from the boundary conditions used in the so-called NASA Langley Zwikker- 
Kosten Transmission Line Code (ZKTL), which finds its origins in the work done by Hersh and Walker 
[2], Heidelberg [3] and Motsinger and Kraft [4]. With slight modifications and rearrangement, the model 
is extended for the present problem to yield the following impedance model in the acoustic domain 

^aBL ~ RqBL + J^aBL ’ (2) 


where the acoustic resistance and reactance contribution from the crossflow boundary layer are, 
respectively 


/? -PS> 

aBL 0 


it— T and X„ L = SSl^- ° 85rf i , 
1.2565/) c 0 C D 1 + 305M.: 


(3) 


where j = V— 1 , p is density, S n is the orifice exit area, C D is the orifice discharge coefficient defined 
below, d is the exit orifice diameter, 8 is the boundary layer thickness, / is frequency, c 0 is the 
isentropic speed of sound, and M ^ is the freestream Mach number. Note that this boundary layer 
impedance model is primarily a function of the Mach number and of the ratio of the orifice diameter to 
the boundary layer thickness. The main contribution of the crossflow is to increase the resistance of the 
orifice and to reduce the effective mass oscillating in the orifice. This model does not provide detailed 
information on how the velocity profile is skewed by the grazing flow. 


The acoustic impedance of the orifice is given by the following expression, where again analytical 
expressions for each lumped parameter are provided in [1], 

^aO — kaolin + ^aOnl + J^^aO * ( 4 ) 


Because of the particular shape of the orifice (aspect ratio h/d <1, beveled orifice), the linear resistance 
R a0 iin ^e to the viscous losses will be small compared to the nonlinear “dump” loss effects that are 
represented by a nonlinear resistance R a0n i . This term is modeled as a generalized Bernoulli flow meter 
and is given by 


R 


pKpQ^ 

2S: 


(5) 


where the dump loss coefficient for the orifice is K D = j ^/l - (/ /c' ; , j , with (> = d/D = the ratio of the 

exit to the entrance orifice diameter, and with the discharge coefficient C D =0.9975- 6.53 (p/Re^) 0 5 , 
Re rf being the Reynolds number based on the orifice exit diameter. For the current configuration, we 
find that K D - 1 . 
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Finally, from Figure 1 the expression for the jet volume flow rate during the expulsion part of the cycle 
can be written in the following form (see Case 1 write-up and [1] for variable definitions): 


Qout Qd 


Z aC + ^aO + ^aBL 


2 C aC {M , 


sAV 

aO + M aB L ) + S ^aC aOlin + R aOnl + ^aBL 


) + l 


( 6 ) 


For an oscillatory channel flow in a circular duct, as discussed in [1], the ratio of the average velocity to 
the centerline velocity is strongly influenced by the Stokes number. Since, for the present case, the 

Stokes number S = ^2nfd 2 /v = 49.47 , the ratio is estimated as u avg /u CL = 0.9 , where u avg corresponds to 
the spatially averaged velocity during the ejection portion of one cycle, and u CL is the centerline velocity. 
Therefore, since LEM gives the volume flow rate through the orifice, the centerline velocity is 


u 


CL 


Qout 

7i(<//2) 2 (0.9) 


(7) 


Implementation and Case Specific Details 

In the approach considered above, the main difficulty comes from the proper modeling of the driver 
volume velocity and cavity volume. Indeed, in this modeling technique the output flow is highly 
dependant on the driver dynamics, especially since at low frequencies, s = jco^ 0 , Eq. (6) shows that the 

response is governed by Q d . It is unclear, however, whether the effective moving area S d of the driver is 
constituted by only the rigid piston or if the flexible membrane displaces some fluid. From a simple 
control volume analysis, assuming incompressible flow and equating the volume fluxes gives 

\Qd\~ 2nfS d W 0 - u CL S n (0.9) = \Q 0Ut \ . (8) 

Substituting in the experimental values that are provided and solving for S d yields an effective area of the 
driver S d = u CLexp S n (0.9)/27i/TF 0 = 1748.5 xlO" 6 mm 2 , which is much less than the value computed from 
the drawings provided on the web site when considering the flexible membrane as a “tensioned drum” 
(s d =4076.7 mm 2 ) . A careful analysis of the piston/membrane dynamics is required to assess the 

severity of this issue. Nevertheless, for the present time, we assume that only the rigid piston displaces 
the fluid, with an area of 2570.4 mm 2 . 


Another point is that the piston-diaphragm is assumed to oscillate in a sinusoidal motion, which turns out 
to be not quite true when comparing the piston displacement over one cycle with a pure sine wave. This 
is actually not unexpected, but it emphasizes the need to model the entire electromagnetic transducer. 


Also, the LEM analysis above can assess the validity of assuming incompressible flow inside the cavity. 
The drive frequency is 150 Hz; it can be shown that a requirement for the incompressibility limit inside 
the cavity is that the operating frequency of the actuator should be much less than the Helmholtz 
frequency. Here the Helmholtz frequency is 


fn 


2 n 


'JC'aC {M a0 + M aBL ) 


- 614Hz . 


(9) 


Next, in order to further characterize the response of the system and to fully understand the effect of the 
nonlinearity present in the orifice resistance, a time domain analysis is performed on the actuator from the 
circuit shown in Figure 1. 


The equation of motion of a fluid particle x out is easily derived and is given by a nonlinear 2 nd -order 
oscillator equation 2 8 3 



1 s 

( M aO +M aBL)^out +R aOnl^ouMout\ + { R aOlin + R uBl) Kut + X out = ~ ~ S ' 11 ) • ( 10 ) 

^ aC ^aC^n 

Similarly, the pressure across the orifice is given by 

A P out =-^W 0 sin ( co?) - p- x out . (11) 

The ODE that describes the motion of the fluid particle at the orifice is numerically integrated using a 4 th 
order Runge-Kutta method with zero initial conditions. The integration is carried out until a steady-state 
is reached. The jet orifice displacement and velocity, pressure drop across the orifice, and the driver 
displacement are shown in Figure 2 for both the (a) linear and the (b) nonlinear solutions of the equation 
of motion. The linear solution is obtained by setting R a0nl = 0 and is performed to verify the physics of 
the device behavior and thus confirm the modeling approach used. The linear solution in Figure 2a shows 
that the pressure inside the cavity (which = the pressure across the orifice A P out ) and the driver motion are 
almost out of phase. All quantities exhibit sinusoidal behavior. The jet orifice velocity x out lags the 
cavity pressure for both the linear and the nonlinear solution. Figure 2b shows the effect of the 
nonlinearity of the orifice resistance. Its main effect is to shift the pressure signal such that the fluid 
particle velocity and the cavity pressure are almost in phase. Also, those two signals exhibit obvious 
nonlinear behavior due to the nonlinear orifice resistance. It is not easy however to compare these plots 
with the experimental data provided. Indeed, we were not clear concerning the relative phasing and data 
processing of the provided experimental data (piston displacement, cavity pressure, velocity 
measurement). 



Figure 2: Time signals of the fluid particle displacement and velocity, pressure across the orifice and 
driver displacement during one cycle: (a) linear solution, and (b) nonlinear solution. The quantities are 
normalized by their respective magnitudes for comparison. 


Unfortunately, our current model does not provide any information on the skewness of the velocity 
profile. Instead, Table 1 below summarizes a comparison between the experiments and the model output. 
Clearly, in its current state, the model overpredicts the output velocity and underpredicts the cavity 
pressure, a result that is entirely possible by overestimating the volume displacement of the driver. Given 
the uncertainty in this quantity, the cavity volume, and the experimental data (not reported) and the 
unknown effect of the small aspect ratio, beveled orifice shape (which violates the model assumptions), 
the discrepancy is not unexpected. We also believe that the boundary layer impedance model is not a 
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significant issue here since at such a low Mach number, it has little effect on the calculated device output. 
Experimental velocity data of the no-flow configuration [M x = 0) is required to address this question. 

Further work is required with regards to the boundary layer impedance model to assess its validity, and it 
must be validated with both numerical and experimental data for a wide range of operating conditions. 
Such work is ongoing. In addition, we are also investigating the effect of the grazing boundary layer on 
the shape or skewness/distortion of the velocity profile. The results of these studies will be reported at a 
later date. 



u c ju x <*> 

A P out (peak-to-peak) [Pa] 

Exp. (LDV) 

E31 

8648 

LEM 

1.72 [Eq. (10)] 

5115 [Eq. (11)] 


Table 1: Comparison between experiments and model 

^ U CL is defined as the maximum amplitude of the phase-locked centerline velocity taken just above the 

orifice center. 
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Introduction 

Separation control by means of steady suction [1] or zero efflux oscillatory jets [2] is known to 
be effective in a wide variety of flows under different flow conditions. Control is effective when 
applied in a nominally two-dimensional manner, for example, at the leading-edge of a wing or at 
the shoulder of a deflected flap. Despite intuitive understanding of the flow, at present there is no 
accepted theoretical model that can adequately explain or describe the observed effects of the 
leading parameters such as reduced suction-rate, or frequency and momentum input. This 
difficulty stems partly from the turbulent nature of the flows combined with superimposed 
coherent structures, which are usually driven by at least one instability mechanism. The ever 
increasing technological importance of these flows has spumed an urgent need to develop 
turbulence models with a predictive capability. Present attempts to develop such models are 
hampered in one way or another by incomplete data sets, uncertain or undocumented inflow and 
boundary conditions, or inadequate flow-field measurements. 

This paper attempts to address these issues by conducting an experimental investigation of a low- 
speed separated flow over a wall-mounted hump model. The model geometry was designed by 
Seifert & Pack, who measured static and dynamic pressures on the model for a wide range of 
Reynolds and Mach numbers and control conditions. [3 ,4] This paper describes the present 
experimental setup, as well as the types and range of data acquired. Sample data is presented and 
future work is discussed. 


Experimental Setup 

The experiment consists of wall-mounted Glauert-Goldschmied type body, [3] mounted between 
two glass endplates where both leading and trailing edges are faired smoothly with a wind tunnel 
splitter-plate (see fig. 1). This is a nominally two-dimensional experiment, although there are 
side-wall effects (3-D flow) near the end-plates. The tunnel dimensions at the test section are 
771mm wide by 508mm high, but the hump model is mounted on a splitter-plate (12.7mm thick), 
yielding a nominal test section height of 382mm (distance from the splitter-plate to the top wall). 
The splitter-plate extends 1935 mm upstream of the model’s leading-edge. The trailing edge of 
the splitter-plate, which is 1 129mm downstream of the model’s leading-edge, is equipped with a 
flap (95mm long), which is deflected 24° upwards to reduce circulation around the splitter-plate 
and avoid separation at the leading-edge. The boundary layer is tripped at splitter-plate leading- 
edge, resulting in a fully-developed turbulent boundary layer (c>~30.5mm) at 2.14 chord lengths 
upstream of the model leading-edge. The tunnel medium is air at sea level. 

The characteristic reference “chord” length of the model is defined here as the length of the 
hump on the wall, i.e. c=420mm. Seifert & Pack [3] used the body virtual leading-edge to define 
their chord length; presently the entire hump length is used as the chord length. As a result of 
this, the current scaled (non-dimensional) coordinates of the overall body shape are slightly 
different from those of [3]. A simple rescaling operation can recover it. 

The model is 584mm wide with endplates at both sides (each approximately 235mm high and 
864mm long). The model is 53.7mm high at its maximum thickness point. Both uncontrolled 
(baseline) and controlled flow scenarios have been considered under various different conditions 
for Re<\, 114,800 and M< 0.12. However, detailed flow field measurements were made at 
7?e=929,000, M=0.100. The boundary layer is subjected to a favorable pressure gradient over the 
front convex portion of the model (fore-body) and separates over a relatively short concave ramp 
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in the aft part of the body. A slot at approximately the 65% chord station on the model, 
immediately upstream of the ramp, extends across the entire span ( s ) of the hump. The model 
was equipped with 165 streamwise and spanwise static pressure ports and 20 dynamic pressure 
ports in the vicinity of the separated flow region. All pressure transducers were calibrated in-situ 
prior to each run. 



Figure 1: Side view of the model and splitter-plate (endplates not shown) 

Separation control is achieved using two methods, namely steady suction and zero efflux 
oscillatory blowing. Both suction and oscillatory blowing are introduced from the spanwise slot. 
Steady suction is achieved by means of a suction pump attached to the plenum with the mass 
flow rate monitored, while zero mass-flux oscillatory suction/blowing is achieved by means of a 
zero efflux actuator specifically designed to minimize three-dimensional effects near the slot. 

Sample Experimental Data 

The primary data acquired for this test case were surface static and dynamic pressures, and two- 
dimensional and stereo (three-dimensional) PIV in the separated and reattachment regions. 
Limited hot-wire and Pitot-tube data was acquired as an independent check of the PIV flow field 
results. An oil-film technique was used to determine the reattachment location. In addition, the 
inflow boundary layer and upper wall (ceiling) boundary layers were documented. 

Baseline Results 


Fig. 2 shows the baseline (no control) surface pressure data, from both dynamic and static 
pressure ports, in the separated flow region. The figure shows that there is no significant 
Reynolds number effect for Re>557,400 on this model. Also the extent of the separated region is 
similar to that of ref. [3] at much higher Re and a different setup and facility (The reference 
pressure in [3] was adjusted by 0.266% in order to match their inflow Cp with the present data). 
The suction peak upstream of the slot, just downstream of x/c= 0.5, is somewhat higher than that 
in ref. [3]. The most probable explanation for this is the difference in the ratio of model height to 
tunnel height for the two cases, namely h/H= 8% [3] versus 13% (present setup). Fluctuation 
pressures showed similar trends for all cases. The flow was shown to be essentially two- 
dimensional in that spanwise pressures did not differ materially in the separated region and 
planes of stereo PIV flow fields in the separated and reattachment regions showed negligible 
spanwise variation (see e.g. fig. 3a and 3b). Oil-film surface shear measurements in the 
reattachment region showed an effectively two-dimensional reattachment line at x/c»l.ll 
(shown on fig. 2). The static and dynamic pressures are virtually insensitive to the presence of 
the slot. This was ascertained by comparing data acquired for the open slot (sealed internally at 
the bottom of the plenum) and the slot sealed externally (not shown). 
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Figure 2. Time-mean and rms surface pressures for the baseline case at various Reynolds 

numbers. 



Figure 3. Spanwise ramp pressures and streamwise velocity from stereo PIV in the vicinity of 
reattachment for the baseline (a,b) and control (c,d) cases. 

Control via Steady Suction 

For the suction test case, control was applied via the two dimensional slot using a suction rate of 
0.01518 kg/s at Re= 929,000. Furthermore, control was applied for the same dimensionless 
conditions at different Reynolds numbers (fig. 4). (Suction rates are often expressed as a mass 
flux coefficient, presently C m =0.15%. Seifert & Pack [3] used C M , to allow direct comparison 
with oscillatory cases.) There is a small Reynolds number effect that can be discerned in fig. 3, 
but the trend is towards the higher Reynolds number data. [3] Additional data acquired at higher 
suction rates (C„~0.456%) showed similar trends to those at lower C /( -0.241%. Near the 
centerline the flow retains its two-dimensional nature (figs. 3c and 3d). 
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PIV Profiles for Baseline & Control Cases 

Examples of two-dimensional (2-D) and stereo (three-dimensional) PIV mean velocity profiles, 
in the vicinity of reattachment, are shown for baseline and steady suction cases in fig. 5. U is the 
streamwise component and V is the component normal to the splitter-plate. Additional 2-D PIV 
data was acquired from upstream of the slot, continuously throughout the reattachment region. 
Based on comparisons with Pitot-tube data, errors associated with 2-D profiles were <1% of the 
maximum velocity, while stereo PIV under-predicted mean velocity profiles by as much as 3% 
of the maximum. 




Figure 5. 2-D and Stereo PIV mean velocity profiles for the baseline (a-d) and control (e-h) cases 
in the vicinity of separation and downstream thereof. 
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Examples of turbulent stresses, corresponding to the velocity profiles above, are shown in fig. 6. 
A preliminary error analysis, based on comparisons with hot-wire anemometry indicates errors 
<20% of the maximum value. 




Figure 6. 2-D and Stereo PIV turbulence profiles for the baseline (a-d) and control (e-h) cases in 

the vicinity of separation and downstream thereof. 


Zero-Efflux Oscillatory Forcing 

A zero-efflux oscillatory jet is produced by a rigid piston, that is secured to the base of the 
plenum by means of a flexible membrane. The piston is driven externally by six voice-coil based 
actuator modules, providing maximum slot velocities of approximately 80m/s at frequencies 
ranging from 60Hz to 500Hz. At the test condition (nominal peak slot velocity = 26.6m/s; 
forcing frequency = 138.5Hz), peak slot velocity measurements vary by less than 3% across the 
span of the slot. Surface pressures and time resolved flow fields are currently being acquired. 
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CASE 3: DETAILS AND SUBMISSION GUIDELINES 


Relevant Details 

The freestream Mach number is M_freestream = 0.10. The atmospheric conditions varied, 
but were essentially standard atmospheric conditions at sea level, in a wind tunnel vented to the 
atmosphere, in a temperature-controlled room. These conditions can be given as approximately: 
p_ambient = approx 101325 kg/(m-s A 2) T_ambient = approx 298 K. Some derived relevant 
conditions are: density_ambient = approx 1.185 kg/m A 3 viscosity_ambient = approx 18.4e-6 
kg/(m-s), u_freestream = approx 34.6 m/s, and Re_freestream = approx 2.23e6 per meter, or 
approx 9.36e5 per chord length of hump 

The upstream boundary conditions from the experiment (associated with the boundary layer 
on the plate at location x/c=-2.14 upstream of the start of the hump, where c = 16.536 inches = 
0.4200 m), to be used to help set/verify CFD inflow BCs, are given in the following files: 

• u-velocity BCs from experiment (pitot probe) 

• u' streamwise turbulent intensity BCs from experiment (hot wire) 

• To complete the Reynolds stress matrix for the upstream BCs at x/c=-2.14, assume 
incompressible fully-developed boundary layer distributions 

There are two required conditions for case 3: (1) no suction case and (2) suction case. A third 
oscillatory control case is also defined; however, there is no experimental data available for it 
yet, so the condition is optional. 

• Condition 1: no flow control (no flow through the 23-inch-span slot; slot left open 
with no suction or blowing through it) 

• Condition 2: suction rate of 0.01518 kg/s (27.13 CFM) through the 23-inch-span 
slot (updated 5 September 2003) 

• Condition 3 (optional): zero-net-mass-flux oscillatory suction/blowing, frequency 
= 138.5 Hz, peak velocity out of slot (Ujetmax) during blowing part of cycle = 26.6 m/s (This 
condition corresponds to an oscillatory blowing momentum coefficient of approximately <cmu> 
= <J>/(cq) = rho*h*(Ujetrms) 2 /(c*0.5*rhoinf*Uinf 2 ) = 0.111%, where h is the slot width) 

Submission Guidelines 

Submissions for this case should include both conditions of no-flow-control and steady- 
suction. The third condition of oscillatory blowing control is optional. The purpose here is to 
determine the state-of-the-art in modeling these types of flows, so we want to explore which 
CFD methods work and which do not. 

• Either model the internal cavity or specify a boundary condition at the slot exit. 

• Either choose to model the flowfield two-dimensionally or three-dimensionally (if three- 
dimensionally, you may choose to model the endplates or you may employ periodic or 
some other appropriate spanwise boundary condition). 

There is the requirement that the details be specific on how the case is modeled, including all 
boundary conditions and approximations made. As the methodologies are assessed at the 
workshop, it will be important to know as many details as possible about the 
calculations/ simulations . 

Detailed requirements include: 

1. The no-flow-control and steady-suction conditions may be run either time- accurately or 
in steady-state mode, depending on your computational method. The optional oscillatory 
blowing control condition (if computed) must be run time-accurately, in order to 
simulate the unsteady nature of the case. 

la. RANS codes run in time-accurate mode (e.g., URANS) solve directly for phase-averaged 

variables, i.e. <Uj> and <u)u )>. Click on link to PDF write-up or link to webpage write-up for 
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more details. Therefore, time-accurate RANS simulations should result in repeating or very- 
nearly-repeating periodic solutions. When periodicity is achieved, averages over one or more 
periods of oscillation yields the long-time-averaged (time-independent mean) values for these 
quantities. 

lb. DNS, LES, or DES simulations will need to be post-processed to obtain the long-time- 
averaged (time-independent mean) values. 

2. GRID STUDY : If the case ismodeled two-dimensionally, then it is strongly suggested that the 
computation be performed using at least two different grid sizes in order to assess the effect of 
spatial discretization error on the solution. If it is modeled three-dimensionally, then solutions 
using more than one grid size are encouraged, but not required. If you use more than one grid, 
submit each set of results separately. 

3. TIME STEP STUDY (for time-accurate computations only): If the case is modeled two- 
dimensionally, then it is strongly suggested that the computation be performed using at least two 
different time step sizes in order to assess the effect of temporal discretization error on the 
solution. If it is modeled three-dimensionally, then solutions using more than one time step are 
encouraged, but not required. If more than one time step is used, submit each set of results 
separately. 

Specific quantities that result from the computations at particular locations will be required for 
submission. Note that for all the following, the coordinate system with x downstream and y up, 
with the (x,y)=(0,0) origin at the start of the hump is adopted. All results from 3-D computations 
are to be given at the z=0 location (centerplane of the tunnel). The chord "c" is the length of the 
bump: 16.536 inches. For the optional condition 3 (oscillatory control), there is only one 
submission requirement: the long-time-averaged Cp (see (e) below). The requirements follow (if 
is not possible to provide a particular quantity, simply leave it out of the "variable" list, and 
reduce the number of columns of data submitted): 

a. x/c vs. Cp on the lower wall for no-flow andsuction cases. For time-accurate computations, 
this should be the long-time-averaged Cp. The definition for Cpis Cp=(p-pinf)/ 
(0.5*rhoinf*uinf A 2), where the "inf" values should correspond to the upstream "freestream" 
location where M=0. 1 . Submit these results for both conditions ofno-flow-through-the-slot, and 
for that of constant suction. Include results at least as far upstream as x/c=-2.14, and at least as 
far downstream as x/c=2.0. Do not include the data on the walls deep inside the slot. 

Name this file: case3.cp.ANYTHING.dat- where "ANYTHING" can be any descriptor you 
choose (should be different for each file if you are submitting multipleruns)-the file should be in 
2-column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #brief description of code/method (pound sign needed) 

6th line: #other info about the case, such as spatial accuracy (pound sign needed) 

7th line: #other info about the case, such as turb model (pound sign needed) 

8th line: #other info about the case (pound sign needed) 

9th line: variables="x/c","Cp" 

10th line: zone t="surface Cp, no flow case" 

subsequent lines: x/c Cp <- this is the data for no flow case 

next line: zone t="surface Cp, suction case" 

subsequent lines: x/c Cp <- this is the data for suction case 

The sample datafile case3.cp.SAMPLE.dat can be downloaded from the website 


b. x/c vs. Cf on the lower wall for no-flow andsuction cases. For time-accurate computations, 
this should be the long-time-averaged Cf. The definition for Cf is Cf=tauw/ (0.5*rhoinf*uinf A 2), 
, where the "inf" values should correspond to the upstream "freestream" location where M=0.1. 
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The term tauw stands for mu*(du/dy) at the wall. Give these results for both conditions of no- 
flow-through-the-slot, and that of constant suction. Include results at least as far upstream as 
x/c=-2.14, and at least as far downstream as x/c=2.0. Do not include the data on the walls deep 
inside the slot. Name this file: case3.cf.ANYTHING.dat -where "ANYTHING" can be any 
descriptor you choose (should be different for each file if you are submitting multiple runs-the 
file should be in 2-column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #brief description of code/method (pound sign needed) 

6th line: #other info about the case, such as spatial accuracy (pound sign needed) 

7th line: #other info about the case, such as turb model (pound sign needed) 

8th line: #other info about the case (pound sign needed) 

9th line: variables="x/c","Cf" 

10th line: zone t="surface Cf, no flow case" 
subsequent lines: x/c Cf <- this is the data for no flow case 
next line: zone t="surface Cf, suction case" 
subsequent lines: x/c Cf <- this is the data for suction case 

The sample datafile case3.cf.SAMPLE.dat can be downloaded from the websit 


c. Velocity profiles and turbulence data at several stations along the lower wall for no-flow and 
suction cases. For time-accurate computations, these should be the long-time-averaged profiles. 
The 15 required stations are: x/c=-2.14, 0,. 2, .4, .65, .66, .8, .9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.6, 2.0. The 

profiles should be taken along a VERTICAL line at each given x-location (not necessarily 
normal to the hump surface), extending from the wall up to at least past the boundary layer 
and/or separation bubble. Also, if your computations included the inside of the slot, include 
profiles along one line at x=0.647 inside the slot (from the lower wall to the upper wall, y goes 
from approximately 0.1 105 to 0.1 142 at this x location). Submit the following quantities: u/Uinf, 
v/Uinf, u'u'bar/Uinf A 2, v'v'bar/Uinf A 2, and u'v'bar/Uinf A 2, where: 
u = horizontal velocity component 
v = vertical velocity component 

u'u'bar = turbulent normal stress in horizontal direction (optional) 
v'v'bar = turbulent normal stress in vertical direction (optional) 
u'v'bar = turbulent shear stress in x-y plane 

Submit two separate files, one for no-flow-control, and one for constant suction. Name these 
files: case3.pro.noflow.ANYTHING.dat case3.pro.suction.ANYTHING.dat-where "ANY- 

THING" can be any descriptor you choose (should be different for each file if you are submitting 
multiple runs)-the file should be in 7-column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #brief description of code/method (pound sign needed) 

6th line: #other info about the case, such as spatial accuracy (pound sign needed) 

7th line: #other info about the case, such as turb model (pound sign needed) 

8th line: #other info about the case (pound sign needed) 

9th line: variables="x/c","y/c", "u/Uinf", "v/Uinf", "uu/Uinf A 2", 

"vv/Uinf A 2","uv/Uinf A 2" 

10th line: zone t="x/c=-2.14" 
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subsequent lines: x/c y/c u/Uinf v/Uin uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=-2.14 

next line: zone t="x/c=0" 

subsequent lines: x/c y/c u/Uinf v/Uin uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=0 

next line: zone t="x/c=0.2" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=0.2 

next line: zone t="x/c=0.4" 

subsequent lines: x/c y/c u/Uinf v/Uin uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=0.4 

next line: zone t="x/c=0.65" 

subsequent lines: x/c y/c u/Uinf v/Uin uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=0.65 

next line: zone t="x/c=0.66" 

subsequent lines: x/c y/c u/Uinf v/Uin uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=0.66 

next line: zone t="x/c=0.8" 

subsequent lines: x/c y/c u/Uinf v/Uin uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=0.8 

next line: zone t="x/c=0.9" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=0.9 

next line: zone t="x/c=1.0" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=1.0 

next line: zone t="x/c=l.l" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=l.l 

next line: zone t="x/c=1.2" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=1.2 

next line: zone t="x/c=1.3" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=1.3 

next line: zone t="x/c=1.4" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=1.4 

next line: zone t="x/c=1.6" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=1.6 

next line: zone t="x/c=2.0" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
at x/c=2.0 

next line: zone t="inside slot, x/c=0.647" 

subsequent lines: x/c y/c u/Uinf v/Uinf uu/Uinf A 2 vv/Uinf A 2 uv/Uinf A 2 <- this is the data 
inside slot at x/c=0.647 

The sample datafile case3.pro.noflow.SAMPLE.dat can be downloaded on the website 
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d. Field line-contour-plots (in one of the following formats: ps, eps, or jpg) of streamlines 
(average streamlines for time-accurate computations) for no-flow and suction cases. These plots 
should be black-and-white line plots only. The plots should go from approximately x/c=0.6 to 
1.6, and should contain enough streamlines to show the approximate size and shape of the 
separation bubble. The x-to-y ratio of the plot should be 1.0. The purpose of submitting these 
plots is to get a qualitative picture of the flowfields. Altogether, submit 2 plots files, one for the 
no-flow-control case and one for the steady suction case. Name these files: 

case3. stream.noflow.ANYTHING.eps 
case3.stream.suction.ANYTHING.eps 

(where the "eps" in this case means encapsulated postscript - use ps, or jpg instead if 
appropriate). 

e. (For optional condition 3 only) x/c vs. long-time-averaged Cp on the lower wall for oscillatory 
case. Include results at least as far upstream as x/c=-2.14, and at least as far downstream as 
x/c=2.0. Do not include the data on the walls deep inside the slot. Name this file: 
case3.cposc.ANYTHING.dat- where "ANYTHING" can be any descriptor you choose (should be 
different for each file if you are submitting multiple runs) -the file should be in 2-column format: 

1st line: #your name (pound sign needed) 

2nd line: #your affiliation (pound sign needed) 

3rd line: #your contact info (pound sign needed) 

4th line: #brief description of grid (pound sign needed) 

5th line: #number of time steps per cycle (pound sign needed) 

6th line: #brief description of code/method (pound sign needed) 

7th line: #other info about the case, such as spatial accuracy (pound sign needed) 

8th line: #other info about the case, such as turb model (pound sign needed) 

9th line: #other info about the case (pound sign needed) 

10th line: variables="x/c","Cp" 

11th line: zone t="surface Cp, oscillatory case" 
subsequent lines: x/c Cp <- this is the data for oscillatory case 

The sample datafile case3.cposc.SAMPLE.dat can be downloaded on the website 
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CASE 3: DIRECT NUMERICAL SIMULATION ON THE CRAY XI 


D. Postl 1 , S. Wernz 1 , and H. F. Fasel 1 

1 Department of Aerospace and Mechanical Engineering, The University of Arizona, Tucson, AZ 85721 


Introduction 

In the present approach, we assess the feasibility of simulating complex flows at high Reynolds numbers 
without the use of turbulence models. With the ever-increasing availability of computational resources on 
supercomputers such as the Cray XI, it has become possible to perform three-dimensional computations 
using hundreds of millions of grid points. Despite the fact that such high Reynolds number simulations may 
not be fully resolved all the way down to the smallest scales, they are considered to be useful for studying 
the development and the dynamics of large coherent structures in complex turbulent flows as well as for 
aiding in our interpretation of results obtained from turbulence models. 


Solution Methodology 

The computer code that was used for the present simulations is based on a highly accurate incompress- 
ible Navier-Stokes solver developed for investigating transitional and turbulent boundary layer flows [1]. 
While the original code was written for a Cartesian grid, it has since been adapted to orthogonal curvilinear 
coordinates [2]. In the code, the 3D Navier-Stokes equations are solved in vorticity-velocity formulation. 
This formulation of the governing equations is obtained by taking the curl of the momentum equation, thus 
eliminating the pressure term. Taking into consideration the fact that both velocity and vorticity vectors are 
solenoidal and using the following definition for the vorticity, 

uj = —V x v, (1) 


the vorticity transport equation 

Du; 1 0 

— - + (v • V) a; — (a; • V) v = v cu + V x F (2) 

ot Re 

is obtained. The velocity field, v, is obtained from the vorticity, cd, through the vector Poisson equation 

V 2 v = V x u. (3) 

The governing equations are non-dimensionalized using a reference velocity Jf e f (e.g. free stream velocity 
Uqq) and a reference length L re f (e.g. cord length c). The volume force, F, on the right-hand- side of 
Eq. (2) will be discussed in the next section. 

For the time-integration of the vorticity-transport equations, a four-stage explicit Runge-Kutta scheme 
with fourth-order accuracy is employed. For the spatial discretization, fourth-order compact differences are 
used in the streamwise and in the wall-normal directions. Assuming periodicity in the spanwise 
direction (z), a pseudo-spectral approach is taken. Each variable is represented by a total of 2K+1 Fourier 
modes: the 2-D spanwise average (zeroth Fourier mode), K symmetric Fourier cosine as well as K anti- 
symmetric Fourier sine modes. As a result of this decomposition, the 3D governing equations reduce to a 
set of 2D equations for each Fourier mode. 
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The equations for the velocity components, Eq. (3), can be reduced to a set of 2 ODEs (for the stream- 
wise and spanwise velocity components) and, as a result of the transformation to orthogonal curvilinear 
coordinates, a 2D convection-diffusion equation with variable coefficients for the wall-normal velocity com- 
ponent (for each Fourier mode). This equation is of the form 

+ « (£, g) • + /? (£, rj) • || + 7 (£, rj) • + S (f , rj) ■ v = f (£, rj) (4) 

and its solution represents the single most computationally expensive step in the numerical integration of the 
governing equations. Eq. (4) is discretized using a fourth-order accurate nine-point stencil. The derivation 
of this stencil is based on an extension to the procedure described in [4]. The resulting linear system is 
solved iteratively using a zebra-alternating-line Gauss-Seidel (ZALGS) algorithm with multigrid accelera- 
tion. Despite its slightly inferior smoothing properties compared to the incomplete-line-LU-decomposition 
algorithm, the ZALGS was chosen because it can easily be vectorized. 

The code has been parallelized using the Message Passing Interface (MPI). The parallelization is done 
in the spanwise direction because the solution of Eq. (4) greatly complicates an efficient implementation 
of typical domain decomposition techniques. Using the present domain decomposition with respect to the 
Fourier modes only yields good parallel efficiency on supercomputers with extremely high communication 
bandwidths. This is due to the fact that the computation of the non-linear terms in the vorticity transport 
equations is done in physical space. The vectorized FFT routines require data to be local to each process, 
which necessitates swapping of the entire 3D arrays before and after each time the flow field is transformed 
from spectral to physical space and vice versa. The Cray XI, with its communication bandwidth of 10 GB/s, 
allows for this task to be performed in a short time compared to the computation of the FFTs. 


Implementation and Case Specific Details 
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Figure 1: Computational grid. Every l(f h point shown in both x and y. 


The grid used for the simulations is shown in Fig. 1. It was obtained from an iterative solver that was 
developed to generate orthogonal curvilinear grids based on the method described in [5]. At each step of the 
iteration, the following set of equations, 


d_ f dxj \ d / 1 dxj \ 

dt, V d£ ) drj \g drj ) 


= 0 , 


i — 1 , 2 , 


(5) 


is solved using a multigrid algorithm similar to the one used in the Navier-Stokes code. The distribution of 
boundary points is obtained from a Newton-Raphson sub-iteration. The only control over grid point cluster- 
ing is given by appropriate specification of the distortion function g. The grid generator was vectorized and 
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Figure 2: Sketch of computational domain. 


domain size (A X/c, A Y/c, AZ/c): 

3.25 , 0.9 , 0.071 

number of grid points in x, y: 

2049, 321 

number of non- symmetric 
Fourier modes in £ (K): 

49 

number of points in £ 

(for computing non-linear terms): 

160 

total number of points: 

~ 105 • W 

Ax + , A y + , A^ + 
(based on x/c = —0.5): 

27... 92, 1.2, 17 

At • Uqo/c (condition 1) 

3.3 • 10- ' 

At • Uqq / c (condition 2) 

1.4 • 10 -4 


Table 1: Computational parameters. 


ported to the Cray XI because, in order to obtain a high-quality orthogonal grid, Eqs. (5) need to be solved 
hundreds of thousands of times. 

An illustration of the computational domain for the Navier-Stokes simulations is given in Fig. 2. At the 
inflow boundary, x/c= —1.05, a laminar Blasius boundary layer profile is imposed with R e% = 22.6 • 10 6 . 
All velocity and vorticity components are prescribed as Dirichlet conditions. For maintaining the fourth- 
order accuracy of the code, the streamwise derivatives are prescribed as well. To prevent reflections from 
the outflow boundary, turbulent fluctuations are damped out in a buffer domain starting at x/c = 1.7 using 
the approach discussed by Meitz and Fasel [1]. At the upper boundary, a slip-wall is imposed, v = 0, and 
irrotational flow is assumed, oj = 0. Since the equations for the streamwise (u) and spanwise (w) velocity 
components reduce to ODEs in the streamwise direction, no upper boundary conditions are required for 
these quantities. At the wall, no-slip and no-penetration conditions are imposed, except over the suction 
slot at x/c « 0.65. The spanwise and streamwise vorticity components at the wall are obtained through 
forward integration by using the known velocity field at the new time level. The wall-normal vorticity 
component is set to zero everywhere except over the suction slot where its value is known from the boundary 
conditions on the velocities. Periodicity is assumed in the spanwise direction with the domain width chosen 
as A Z/c = 0.071. This seemingly narrow domain was chosen as a compromise to achieve a reasonable 
spanwise resolution. The choice was also based on the authors’ conviction that the coherent structures 
present in this flow are predominantly two-dimensional in nature. The relevant computational parameters 
are summarized in table 1 . 

The laminar boundary layer is tripped to turbulence near the inflow boundary (see Fig. 2) by introducing 
high-amplitude, 3D disturbances into the flow using a volume forcing technique, i.e. by adding a time- 
harmonic forcing term, V x F v f, to the right-hand- side of the vorticity transport equations (2). The force is 
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applied for selected spanwise Fourier components, k. For each of these k , it takes on the form of 


N 

F$(£,77,M) = ^7), (6) 

n = 1 

where vector A^ k ^ = (A^ A V ,A Z ) determines the amplitude and spatial orientation of the force, and / n,A 0 

and 0^^ are the frequency and phase angle in time. The spatial distribution of the volume force, ^ k \ 
has the Gaussian shape 


q(n,k) 

*vf 


exp 


t _ An,k) y 

S W 

(n,k) 

a vf 


(n,k ) ' 


An,k) 

V 


(7) 


with and b^ k ^ determining its size. The wall-normal forcing locations were chosen such as to achieve 
maximum receptivity of the flow to the disturbance input. Due to the forcing, the flow transitions to turbu- 
lence and at approximately x/c = —0.5, a fully developed turbulent profile is reached. 

As a result of the tripping procedure described above, it was found to be difficult to precisely match 
a particular turbulent profile at a given streamwise location. Due to the limitations on available compu- 
tational resources, it was impossible to trip a laminar boundary layer, match the given turbulent profile at 
x/c = —2.14 and continue the computation all the way to and beyond the hump. Since it was noted by 
Seifert and Pack [6] that, in their experiments, the thickness of the upstream boundary layer had only a 
minor effect on the flow field, it was deemed more important to reach a fully turbulent state upstream of 
the hump than to precisely match the experimental boundary layer. In the present simulations, the approach 
boundary layer is thinner than the one in the experiments. 

For the case involving steady suction (condition 2), the slot was modeled by a boundary condition on 
the wall-normal velocity component, 


Vsuct — A t 


suet 


cos ( 2n 


Wslot / J 


1 3 


( 8 ) 


where s represents the distance from the center of the slot along the body surface. The width of the slot, 
w s iot , was increased to obtain a reasonable spatial resolution (see Fig. 2). The strength of the suction, 
A suc t = —0.3, was adjusted accordingly to match the specified suction rate of 0.015 kg/s. 




Figure 3: Instantaneous gray-scale contours of the spanwise vorticity. Unforced case. 
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Instantaneous gray-scale contours of the spanwise vorticity component are shown in Fig. 3 for the un- 
forced case (condition 1) and in Fig. 4 for the case with steady suction (condition 2). 
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Introduction 

The development of computational tools that can be used to guide technologies for flow control applications 
at realistic Reynolds numbers and in complicated configurations is a topic of significant interest. Compu- 
tational Fluid Dynamics (CFD) offers a useful tool in understanding flow characteristics and studying the 
performance gains that can be achieved through flow control, though advances in several areas are needed 
to improve the robustness of CFD predictions for applications. 

Most of the flow fields for which control would be advantageous are often very complex - the flows 
are far from equilibrium, e.g., separated or on the verge of separation, and distorted by effects such strong 
pressure gradients and streamline curvature. Application Reynolds numbers are sufficiently high that em- 
pirical input to any simulation strategy appears unavoidable. In addition, it would be desirable in many flow 
control applications to exert a “microscopic” (e.g., small-scale) input and observe a desired “macroscopic” 
(large-scale) output. This implies a wide range in the geometric scales to be simulated, which in turn places 
significant burdens on issues such as grid design and construction, possibly more so than in most conven- 
tional CFD applications. Thus, the opportunity for a thorough assessment of simulation strategies against 
measurements and other simulation approaches represents a valuable undertaking. Summarized in this re- 
port are predictions of Case 3 using solutions of the Reynolds-averaged Navier-Stokes (RANS) equations 
and Detached-Eddy Simulation (DES). 



Figure 1: Contours of the instantaneous vorticity magnitude in the leeward region of hump. Flow field 
prediction of the no-flow-control case obtained using DES. Surfaces colored by pressure. 
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Solution Methodology 


The compressible Navier-Stokes equations are solved using Cobalt. The numerical method is a cell-centered 
finite volume approach applicable to arbitrary cell topologies (e.g, hexahedra, prisms, tetrahedra) and de- 
scribed in Strang et al.[ 1]. The spatial operator uses the exact Riemann Solver of Gottlieb and Groth [2], 
least-squares gradient calculations using QR factorization to provide second order accuracy in space, and 
TVD flux limiters to limit extremes at cell faces. A point implicit method using analytic first-order inviscid 
and viscous Jacobians is used for advancement of the discretized system. For time-accurate computations, a 
Newton sub-iteration scheme is employed, the method is second order accurate in time. The domain decom- 
position library ParMETIS [3] is used for parallel implementation and provides optimal load balancing with 
a minimal surface interface between zones. Communication between processors is achieved using Message 
Passing Interface. 


Model Description 

Summarized in the next section are the simulation parameters for both two-dimensional (2D) and three- 
dimensional (3D) configurations using RANS and DES. RANS predictions are obtained using two turbu- 
lence models: the Spalart-Allmaras one-equation model [4] (referred to as S-A throughout) and the two- 
equation SST model of Menter [5]. All of the calculations summarized below are of fully turbulent flows, 
i.e., with turbulent boundary layers initiated along all solid surfaces of the computational domain. 

The DES formulation is obtained by replacing in the S-A model the distance to the nearest wall, d, by 
d , where d is defined as, _ 

d = min(d, Cdes A ) , Cdes = 0.65 . (1) 

In Eq. (1), A is the largest distance between the cell center under consideration and the cell center of the 
neighbors (i.e., those cells sharing a face with the cell in question). The location whered is determined by 
the grid spacing, i.e., d = Cdes A, defines the interface between the RANS region and the LES region. In 
applications for which the wall-parallel grid spacings (e.g., streamwise and spanwise) are on the order of 
the boundary layer thickness, the RANS region comprises most or all of the boundary layer and the closure 
applied is the S-A RANS model. In the LES region the closure is a one-equation model for the subgrid scale 
eddy viscosity. For the present configuration, the combination of the adverse pressure gradient upstream 
of the hump thickened the boundary layer. That feature combined with the spacings for the current grids 
caused the RANS-LES interface to reside sufficiently within the boundary layer such that the upstream 
flow prediction was less accurate compared to RANS results and experimental measurements. To provide 
an evaluation of the DES predictions against the RANS for nominally similar upstream conditions, RANS 
behavior was maintain to x/C = 0.65, slightly upstream of the slot (C is the hump chord length). 


Implementation and Case Specific Details 

All of the computations to be reported at the workshop are summarized in Tables 1-3. The table reports the 
grid sizes, grid topology, whether the computation was 2D or 3D, the boundary conditions, and turbulence 
models employed for a given simulation (the “X” in the table indicating the simulation was performed). 
The nomenclature in the tables for the grid topology indicates whether the mesh was structured or unstruc- 
tured and with additional details for the 3D computations as summarized below. The nomenclature for 
the boundary conditions “slip top wall” indicates a slip condition was applied to the upper surface of the 
computational domain, “bl top wall” indicates a boundary layer grid on the top wall and that the no-slip 
condition was applied, and “bl all walls” indicates boundary layer grids on all walls and the imposition of 
no-slip conditions on all solid surfaces. The distance from the lower horizontal surfaces (that are faired to 
the upstream and downstream sides of the hump) to the upper boundary of the computational domain was 
15.032 inches (slightly less than one chord length). 
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grid size 

grid 

topology 

boundary 

conditions 

RA 

S-A 

jsts 

SST 

DES 

S-A 

2D 

421 x 109 

structured 

slip top wall 

X 

X 


841 x 101 

structured 

slip top wall 

X 

X 


841 x 217 

structured 

slip top wall 

X 

X 


841 x 257 

structured 

slip top wall 

X 



841 x 257 

structured 

bl top wall 

X 



890 x 257 

structured 

bl top wall 

X 



1.14 x 10° 

unstructured 

slip top wall 

X 



2.47 x 10° 

unstructured 

slip top wall 

X 



3D 

841 x xlOl x 41 

structured, periodic 

slip top wall 



X 

2.59 x 10° 

unstructured, half-geometry 

bl all walls 

X 



4.90 x 10 b ' 

unstructured, half-geometry 

bl all walls 

X 

X 


10.72 x 10 b 

unstructured, half-geometry 

bl all walls 

X 



9.14 x 10° 

unstructured, full-geometry 

bl all walls 



X 


Table 1 : Summary of the computations of the no-flow-control case. 


Parameters summarized in the tables correspond to three sets of calculations. The predictions of the “no- 
flow-control” case in Table 1 were performed using configurations that meshed the cavity and slot, though 
with the lower (horizontal) cavity surface closed. The steady- suction case imposed a fixed mass flux of 

0.01518 kg/s (divided by the 23 inch span of the slot through which suction was applied in the experiments) 
along the entire lower cavity surface. The zero-net-mass-flux oscillatory suction/blowing was performed 
by prescribing a sinusoidal variation in the mass flux exiting/entering the lower cavity surface. The suc- 
tion/blowing frequency was that prescribed on the workshop web page, i.e., 138.5 Hz. The peak mass flux 
was adjusted at the lower cavity surface in order that the maximum velocity out of the slot was approxi- 
mately 26.6 m/s, which corresponded to a mass flux of 0.0179 kg/s (over the entire slot) as summarized 
in Table 3. As also shown in Table 3, additional suction/blowing rates were applied in order to investigate 
the characteristics of the jet velocity through the slot and impedance effects associated with large blowing 
coefficients. The results from these simulations will be presented at the workshop, though are not included 
with the data files that were submitted. 

For each simulation, the reference conditions corresponded to standard atmosphere. At the outlet of 
the computational domain the pressure was prescribed at 14.696 psia, the reference temperature was 519 
Rankine, and the corresponding density specified using the ideal gas law. The reference Mach number was 
0.1, leading to a Reynolds number per unit of length of 5.899 x 10 4 , which yields a chord-based Reynolds 
number of 9.75 x 10 5 for the simulations summarized below. For the 2D simulations the reference pressure 
used in calculation of the pressure coefficient was adjusted in order to match experimental measurements 
upstream of the hump leading edge. For the 3D calculations the pressure coefficient was computed using 
the reference conditions in the freestream at x/C = —2.14 (the coordinate origin x/C = 0 corresponding 
to the hump leading edge). 

Additional notes on the 2D cases: 

1. RANS calculations were performed on two of the 2D structured grids available on the workshop web 
site and on two of the 2D unstructured grids available on the workshop web site. The 2D structured 
grids comprised of 841 x 217 cells and 421 x 109 cells are referred to as “Structured 2D Grid #1” and 
“Structured 2D Grid #2” on the web site. The unstructured grids comprised of 2.45 x 1(? cells and 
1.12 x 10 5 cells are referred to as “Unstructured 2D Grid #1” and “Unstructured 2D Grid #2” on the 
web site. 

2. The structured grid comprised of 841 x 101 cells was created using Gridgen, and starting with the 
“Structured 2D Grid #1”. The grid resolution was coarsened in the wall normal direction, while the 
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grid size 

grid 

topology 

boundary 

conditions 

RA 

S-A 

lNS 

SST 

DES 

S-A 

2D 

421 x 109 

structured 

slip top wall 

X 

X 


841 x 101 

structured 

slip top wall 

X 

X 


841 x 217 

structured 

slip top wall 

X 

X 


3D 

841 x 101 x 41 

structured, periodic 

slip top wall 



X 

4.90 x 10 b 

unstructured, half-geometry 

bl all walls 

X 

X 



Table 2: Summary of the computations of the steady- suction case. 



grid size 

grid 

topology 

boundary 

conditions 

RANS 

S-A 

mass 

flux (kg/s) 

2D 

841 x 101 

structured 

slip top wall 

X 

0.0179 

841 x 101 

structured 

slip top wall 

X 

0.0140 

841 x 101 

structured 

slip top wall 

X 

0.0210 

841 x 101 

structured 

slip top wall 

X 

0.0280 

841 x 101 

structured 

slip top wall 

X 

0.0560 


Table 3: Summary of the computations of the pulsed blowing and suction case. 


surface point distribution was left unchanged. The first ^/-location off the wall was also unchanged. 
This grid was created in order to provide a more efficient mesh for subsequent DES calculations that 
extruded this grid into the span. 

3. With the exception of the 2D RANS calculation of the baseline case using the grid comprised of 
890 x 257 cells, the upstream section of the computational domain extended 6.39 chords forward of 
the hump leading edge. The simulation domain for the case with the grid size 890 x 257 extended 
10.64 chord lengths upstream of the leading edge. The purpose of the computation was to investigate 
the influence of thicker upper-wall boundary layer on the pressure distribution over the lower surface 
(on which the hump is mounted). Along the lower surface of the domain from 10.64 to 6.39 chord 
lengths a slip condition was applied, with a no-slip condition applied at the streamwise location 6.39 
chord lengths upstream of the hump leading edge. This in turn allowed the boundary layer to develop 
from the same upstream location as the other simulations. 

4. The downstream section of the domain for all 2D computations extended three chord lengths from the 
trailing edge of the hump. 

5. With the exception of the simulations performed of the case with pulsed suction/blowing (c.f., Ta- 
ble 3), all RANS predictions are of the steady state flow. The governing equations were integrated 
using large timesteps, corresponding to a CFL = 10P. 

6. The predictions of the pulsed suction/blowing case were time-accurate. Based on the results of a 
timestep study, the timestep non-dimensionalized by the hump chord length and freestream velocity 
(computed using the reference Mach number of 0. 1 and reference temperature) was 8 x 1CT 4 (1 x 10 -5 
seconds). The mass flux values in Table 3 correspond to the peak values (divided by the 23 inch span 
of the slot). 

Additional notes on the 3D cases: 

1. All DES predictions were of the unsteady flow using time-accurate computations. A timestep study 
was performed using the DES of the baseline configuration (no flow control) and the grid comprised 


3.4.4 





of 841 x 101 x 41 cells. Dimensional timesteps of 2 x 10 -5 , 4 x 10 -5 , and 8 x 10 -5 seconds (corre- 
sponding to dimensionless values of 0.0016, 0.0032, and 0.0064) were employed in the calculations. 
Evaluation of frequency spectra at three locations in the separated region indicated relatively small 
differences in the spectra for simulations using timesteps of 2 x 1CT 5 and 4 x 10 -5 seconds. Con- 
sequently, for the DES predictions summarized in Tables 1 and Table 2 the timestep was 4 x 1CT 5 
seconds. 

2. Tests of the time-dependent nature of the RANS flow fields were performed for a few select cases 
in which the RANS equations were integrated in a time-accurate fashion. These tests showed that 
the time-accurate RANS predictions evolved to steady solutions. Thus, the 3D RANS predictions 
summarized in the tables and to be presented at the workshop are of the steady state flow, obtained by 
integrating the equations using large timesteps, corresponding to a CFL = ltf . 

3. The DES prediction in Table 1 obtained using the structured grid was of a section of the hump. A 
grid with 841 x 101 cells (and used in the RANS calculations as summarized in Tables 2 and 3) was 
extruded into the spanwise direction to create the 3D geometry. The spanwise dimension was meshed 
using 41 points with an equal spacing of 0.05 inches, leading to a spanwise period of 0.12 chord 
lengths. Periodic boundary conditions were applied along the spanwise direction, indicated by the 
reference “structured, periodic” in the grid topology entry in the tables. 

4. Aside from the structured grid, the 3D computations summarized in Table 1 and Table 2 considered 
the entire tunnel geometry, i.e., from the splitter plate to the upper wall and spanwise extent of the 
lower test surface, including the endplates that are attached to the hump. The RANS predictions are of 
half the geometry with a symmetry condition imposed along the centerplane (indicated by the “half- 
geometry” entry in the tables). The DES prediction on the grid with 9.14 x 1CP cells is of the entire 
geometry, i.e., without any symmetry conditions imposed (indicated by the “full-geometry” entry in 
the tables). Contours of the instantaneous vorticity from the DES prediction are shown in Figure 1. 
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Introduction 

The primary purpose of this study is to compare a high-order compact scheme with standard second-order 
methodologies for flow over a hump. The secondary purpose is to explore predictive capabilities of high- 
order hybrid techniques with those of standard Reynolds-Averaged Navier-Stokes (RANS) approaches. 


Solution Methodology 

The governing equations are the unsteady, three-dimensional, full Navier-Stokes equations. The governing 
equations are expressed in general curvilinear coordinates and cast into strong conservation-law form. All 
dependent variables have been normalized by reference values except for p, which has been nondimensional- 
ized by p re fU 2 e j. Sutherland’s law for the molecular viscosity coefficient p and the perfect gas relationship 
are also employed. Stokes’ hypothesis has been invoked for the bulk viscosity coefficient. 

Two forms of the above equations are implemented into the current scheme, the unfiltered and the 
Reynolds-averaged forms of the Navier-Stokes equations. The unfiltered form of the equations are em- 
ployed for laminar, large eddy (LES) and direct numerical simulations (DNS). Standard compressible LES 
approaches use Favre-filtered Navier-Stokes equations, which result in additional subgrid-scale stress and 
heat flux terms that must be modeled with subgrid- scale (SGS) models. The present approach, referred to 
as the implicit LES (ILES) methodology, uses the unfiltered form the the equations and the scheme’s filter in 
lieu of a SGS model to prevent pile up of energy at the high- wave numbers of the mesh. Classical Reynolds 
averaged simulations (RANS) are used to account for turbulence effects through the incorporation of an 
eddy viscosity model. The two-equation k — e model has been employed to compute the eddy viscosity. 
To better improve flow simulations in separated flow regions the k — e model has been modified to blend 
between the RANS model and the ILES approach. Hereafter this is referred to as the high-order hybrid 
approach. 

Time Integration 

The solver advances the solution with second-order temporal accuracy using an implicit, subiterative Beam- 
Warming algorithm [1]. The approximate-factorization form of the Beam- Warming algorithm [1] may be 
written in delta form as 
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where d[F]/dQ, d[G\/dQ , and d[H]/dQ, are flux Jacobians, Q is the solution vector, and S represents the 
spatial difference operator. Newton-like subiterations are incorporated to reduce linearization and factoriza- 
tion errors thereby improving the temporal accuracy and stability properties of the algorithm. Subiterations 
also permit the use of the more efficient diagonal form [2] of the implicit algorithm while retaining time 
accuracy. The subiterations are given by A Q = QP +1 — Q p where p denotes the number of subiterations. 
For the first solver application (p = 1), QP = Q n , where n is the solution time-level. With subsequent subit- 
erations Q p+1 — >► Q n+1 . Based on previous unsteady flow computations, three applications of the solver per 
time step are typically found to be sufficient in order to recover second-order time accuracy. 

The implicit portion of the algorithm uses second-order accurate three-point backward differencing for 
the time derivative and second-order central finite-difference approximations for the spatial derivatives (de- 
noted with the subscript (2) in Eqn. 1. Nonlinear artificial dissipation (not shown) is appended to the 
implicit operator to enhance stability. Note that the spatial derivatives on the right-hand side of Eqn. 1 use 
high-order compact-difference operators (described below) that recover higher-order spatial accuracy with 
subiterations. 

The k — e equations are solved uncoupled from the flow equations with a procedure similar to the one 
described above. Since the flow equations are solved uncoupled from the turbulence model, subiterations 
not only have the benefits described above, but also eliminate the time lag between the flow and turbulence 
model equations. 


Spatial Discretization 

The governing equation spatial derivatives are discretized using the compact finite-difference scheme of 
Lele [3]. Using the tridiagonal subset of this scheme, the spatial derivative of any scalar /, such as a metric, 
flux component, or flow variable, can be obtained by solving the system: 


a 




T - ( y. 


df_\ _ ( fi+l - fi—1 

dUi ° 


i + 1 


L) + » ( *+»-*-» ) 


( 2 ) 


where a , a and b determine the spatial properties of the algorithm. Of all the schemes this system generated, 
only the the compact fourth-order scheme (C4), corresponding to a a = |,6 = 0, was used in the 
present investigation. 

Dispersion-error characteristics and truncation error of these schemes are discussed in detail in Refs. 
[3] and [4]. Also discussed in Ref. [4] are descriptions of the higher-order one-sided formulas used near 
boundaries. 


Low-Pass Filter Scheme 

The non-dissipative compact difference scheme incorporated a non-dispersive spatial filter to eliminate nu- 
merical instabilities [5]. The filter is applied to the solution vector Q in each of the three computational 
directions following each subiteration. The filtered values of the solution vector are obtained by solving the 
tridiagonal system 


OLfQi - 1 + Qi + OtfQi+l — ^ — (Qi+n + Qi-n ) ( 3 ) 

n = 0 Z 

where Q represents the filtered value of Q. aj is a free parameter for this family of filters which must remain 
in the range —0.5 < aj < 0.5. Higher values of aj correspond to less dissipative filters. The coefficients 
a n in Eqn. 3 are summarized in Ref. [4]. On uniform meshes, these filters preserve constant functions, 
completely eliminate the odd-even mode decoupling, and do not amplify waves [6]. Low-pass filtering 
provides dissipation at the high modified wave numbers only where the spatial discretization already exhibits 
significant dispersion errors, whereas non-spectral based artificial viscosity and upwind-biased schemes 
introduce dissipation across a wide range of wavenumber. 
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Filtering of the first four points near any domain boundary requires the implementation of one-sided 
high-order Pade-type formulas as discussed in Refs. [4] and [5]. Unlike the interior filter, the one-sided high- 
order filters can amplify wave amplitudes for certain ranges of wavenumber. To eliminate this undesirable 
amplification behavior, higher values of aj are used near boundaries [8]. 

Parallelization and Multidomain Interface Treatment 

The above components are implemented in a parallel scheme based on an overset grid methodology. The 
overset grid implementation uses sixth-order explicit Lagrangian interpolation for the overlapping blocks. 
The overset grid methodology is used in conjunction with the MPI library for inter-block communication. 
The resulting MPI code has been successfully ported to several parallel platforms. Additional background 
on the parallel compact solver is covered in detail in Ref. [9]. 


Model Description 

The k — e model has been chosen to calculate the eddy viscosity, , for RANS and high-order hybrid 
RANS/LES solutions. The Standard Jones and Launder model with the Low-Reynolds-Number corrections 
of Launder and Sharma [10] is employed as described in Ref. [11]. The k — e equations, like the flow 
equations, are also expressed in general curvilinear coordinates and conservation-law form. The current 
implementation varies from the standard model in the production term Fj^ which is now limited to prevent 
anomalous behavior in regions where the the eddy viscosity should be near zero [11]. The term i/J is now 
given by 

Pk = min ( P k , 20. Ope) (4) 

where 

p k = (il) (</>- I (lit) 2 ) ^ = M + 2min(0,|a| - | w |) (5) 

1 5 1 is the magnitude of the rate of strain tensor and \cu\ is the magnitude of the vorticity. This formulation 
reduces the eddy viscosity in regions of the flow where vorticity exceeds the rate of strain, such as in a vortex 
core, and has minimal effect in shear layers. 

The high-order hybrid RANS/LES approach consists of two modifications to the k — e model described 
above. The first modification defines a new dissipation length scale for the k — e model based on the 

minimum grid cell size. The quantity [—pe] in the source term of the k equation is replaced with [ pk / ] 
where 


l = min (7/c-e, CA ) , l^_ e = C = 0.61, A = max (A s x , As y , A s z ) (6) 

The second modification applies a double blending function which turns off the eddy viscosity away 
from the wall and in the separated flow region downstream of the hump. The blending function in the wall 
normal direction is given by 


d 

aX 

where pt-LES is zero for the ILES approach, A = fc 3 / 2 /e, d is the distance to the nearest wall, and a is a 
constant chosen such that T = 0.5 in the log-law region for a flat plate boundary layer. 

The streamwise blending function uses a sum of two hyperbolic tangent functions which rapidly turn off 
the eddy viscosity in the region between the downstream lip of the slot and x/c = 1.5. The total blending 
function is a product of the streamwise and normal blending functions. 



Pt = (1 - r)p t -RANS + r pt-LES, r = tanh f 
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Table 1: Case Summary 


Case 

Grid 

Dimension 

Scheme 

Order 

Model 

At 

Span 

z/c 

Flow Control 

None 

Suction 

Oscillatory 

i 

Coarse 

2D 

2nd 

RANS 

i.o x Hr 4 

n/a 

X 

X 

X 

2 

Fine 

2D 

2nd 

RANS 

5.0 x 10“ 5 

n/a 

X 

X 

X 

3 

Coarse 

2D 

4th 

RANS 

1.0 x 10 -4 

n/a 

X 

X 

X 

4 

Coarse 

3D 

2nd 

RANS 

1.0 x 10 -4 

0.36 

X 

X 


5 

Coarse 

3D 

4th 

RANS 

1.0 x 10 -4 

0.36 

X 

X 


6 

Coarse 

3D 

4th 

HYBRID 

1.0 x 10 -4 

0.13 

X 

X 



Implementation and Case Specific Details 

Multiple simulations were completed for this workshop which are summarized in Table 1. This study 
employed modified versions of the fine grid provided on the workshop website. The four-block mesh was 
reduced to a two-block mesh by eliminating the fourth block upstream of x/c = —2.14 and by joining the 
second and third blocks which make up the slot and cavity used for flow control. The two-block grid was 
then smoothed while enforcing orthogonality at the walls with an elliptic solver. Finally, the slot mesh was 
extended five points into the hump grid to produce an overset topology. A coarse mesh was derived from 
this modified fine grid by eliminating every other point in each coordinate direction. For two-dimensional 
parallel calculations, the coarse mesh was divided into 19 blocks and the fine grid into 36 blocks. 

Two three-dimensional meshes were also developed from the coarse 2 D mesh. Both employed 51 points 
in the spanwise direction. The first 3 D grid assumed symmetry at the span centerline so the half-span wind 
tunnel mesh had a A z/c = 0.0139. The second mesh, used for periodic computations, employed a smaller 
span of approximately two boundary layer heights at the inflow corresponding to 0.13304c which led to a 
A z/c = 0.00289, when the periodic five-point overlap was removed. The 3 D grid was partitioned into 56 
blocks for the simulation. 

Boundary conditions used in the simulations include adiabatic, no-slip walls for the lower surface of the 
wind tunnel, the hump, and walls of the slot. Inviscid wall boundary conditions were imposed for the top of 
the wind tunnel and the bottom of the cavity. For the no-suction case, v is set to zero at the bottom of the 
cavity. Similarly, suction and oscillatory blowing and suction cases are modeled by prescribing v at the base 
of the cavity to match the appropriate mass flow rate or slot velocity. Velocities and density are prescribed at 
the tunnel inflow, and pressure is updated with a Neumann boundary condition (i.e. ^ = 0 ). This inflow 
profile, generated from a separate flat plate RANS calculation, matches the momentum thickness computed 
from the given experimental inflow profile. At the outflow, pressure is fixed and all other variables are 
updated using Neumann boundary conditions. For the two three-dimensional computations, symmetry and 
an inviscid wall boundary conditions are applied at the spanwise boundaries for the cases that model half 
the wind tunnel span, and periodicity is imposed in the spanwise direction of the second thin 3 D mesh. 

Due to the multiple parameters and cases being simulated, a single approach was not used to initialize 
the flow. In general, the two-dimensional coarse mesh cases initialized the hump domain with the inflow 
profile and specified no flow in the slot and cavity. The cavity pressure, density, and turbulence variables (k, 
and e), were set to the wall values of the inflow profile. Since the coarse and fine grid locations are identical 
at every other grid point, the two-dimensional fine grid flow variables were initialized using the coarse grid 
solution at coincident grid points. The non-coincident mesh points were initialized by linearly interpolating 
the flow variables of surrounding points. Three-dimensional cases were initialized using the corresponding 
2D solution for each spanwise plane. 

The major findings of the study are 

• The pressure coefficient was qualitatively correct for all cases. The peak magnitude was slightly 
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smaller than experimental values downstream of the hump mid-span and the separated flow region. 


• Three-dimensional RANS simulations of the entire wind tunnel did not develop any appreciable span- 
wise flow variation. 

• Preliminary computations with a Hybrid approach performed worse than the traditional RANS model. 
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Introduction 

The Flow Simulation Methodology (FSM) is applied to the proposed experimental geometry. In this ap- 
proach only the largest scales of motion need to be resolved, while in regions where there are no large 
coherent turbulent structures a state of the art RANS model is recovered. This allows us to obtain good 
results with fewer than three million grid points for a 3-D calculation. 

The simulations performed for this workshop build on our extensive previous work with this geometry, as 
well as our simulations of the experiments of Seifert and Pack [1, 2, 3, 4]. We have performed both laminar 
[5] and turbulent [6] simulations for the Seifert and Pack geometry. Lower Reynolds number (Re = 1(f) 
turbulent cases have been computed in both 2-D and 3-D. Additionally, 2-D simulations were performed at 
Re = 10 6 . The same geometry was also used to test closed loop control strategies [7]. For these earlier 
investigations a high order (fourth-order in time and space) explicit finite difference code was employed. 
This code was a curvilinear coordinate extension of previous codes developed for transition research, and 
consequently was designed to accurately predict the growth rates of disturbances. 

For this CFD challenge the computational cost of using an explicit code, especially in 3-D, was unac- 
ceptably high. For this reason we have implemented our FSM approach [8, 9] in the CFL3D code developed 
at NASA Langley. This is a second order implicit thin-layer Navier-Stokes code. All the results submitted 
for this workshop were performed using this code, with the viscous terms activated in all three coordinate 
directions. 


Solution Methodology 

Details regarding CFL3D can be found in the CFL3D User’s Manual [10]. Implementation of the FSM 
(described below) does not require changes to the solution methodology of CFL3D. Since the FSM is a time 
dependent model, all simulations were run in an unsteady fashion using the second order r — TS pseudo-time 
iteration scheme. 


Model Description 


In order to bridge the gap between RANS and LES Speziale [11] proposed a methodology for computing 
turbulent flows, which we later called Flow Simulation Methodology. In the FSM, in which a contribution 
function is introduced which compares the local grid resolution relative to some local turbulent length scale 
in order to provide the correct magnitude of the turbulent stresses. The contribution function is designed 
so that in the limit of high resolution the Reynolds stresses vanish, yielding a DNS. For the limit of low 
resolution, the contribution function goes to unity, yielding a RANS. The subgrid stress is then computed as 


T ij — / A R 


RANS 

ij 5 


( 1 ) 


where /a is the contribution function. 


3.6.1 




Figure 1: Instantaneous iso-surface of total vorticity. 


The FSM differs from typical hybrid or zonal models in two important ways. First, since the contribution 
function is a time and space varying function, there is no need to specify a priori which regions of the flow 
are computed as RANS and which as LES or DNS. Second, since a single underlying RANS model is used, 
the FSM is consistent in its limits; as opposed to zonal models which involve the ad-hoc patching of two 
unrelated models. 

For the simulations presented here, the underlying RANS model is an Explicit Algebraic Stress Model 
(EASM) [12, 13]. To obtain the required turbulent length scales, the k — e equations are solved. No ramping 
functions are used in the evaluation of the eddy-viscosity. However, in order to eliminate the near wall 
singularity in the £ equation an f £ 2 damping function is employed [14]. At present, the FSM implementation 
in CFL3D is only applied to the linear eddy-viscosity terms so the non-linear terms in the EASM are not 
used in these simulations. 

Several variant forms of the contribution function have been proposed. In the present case the following 
form is used: 


f(A,L,L K ) 


CpA 2 / 3 - L ^ /3 
L 2 /3 _ L 2 ^ 


where A = ( AxAyAz ) 1 / 3 , and the turbulent length scales were estimated using L = Lk = 

(/i 3 /^) 1 / 4 . Cd is a constant which roughly represents the number of grid points per resolved structure. 
Clipping is applied to insure that 0 < /a < 1. This function is designed so that the RANS limit is reached 
as the grid spacing approaches the integral scale, and the DNS limit as the grid spacing approaches the 
Kolmogorov length scale. The specific form is motivated by a Kolmogorov inertial scale five-thirds decay, 
however this should not be considered a rigorous justification of this form since we do not expect to see a 
significant inertial spectrum in the flows we are investigating. 

A typical instantaneous flow field is shown in figure 1. One can clearly see the development of large 
turbulent structures in the recirculation region whereas the oncoming turbulent boundary layer is steady. 


Implementation and Case Specific Details 

A nearly orthogonal grid is employed which consists of a two-dimensional body fitted grid over the model 
which is then extruded in the spanwise dimension. The grid is generated using an elliptic equation [15]. 
The method is equivalent to the composition of a quasi-conformal mapping and a non-uniform rectilinear 
stretching. Grid point clustering is achieved by modifying the two independent stretching functions (one in 
each grid index). The stretching is chosen to cluster points near the wall, and to increase the resolution over 
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Figure 2: Computational grid. Every 4th point shown in x and y. 


the concave portion of the ramp. (The quasi-conformal map alone would yield a grid with few points in 
concave regions.) The grid is shown in figure 2. 

The computational domain consists of a portion of the test section extending one chord length upstream 
of the model, and two chord lengths downstream. The lower surface is a no-slip constant temperature wall. 
In order to avoid the need for additional grid points, the boundary layer on the upper wall is not simulated. 
Instead an inviscid surface is used. 

At the inflow a turbulent boundary layer is prescribed. For computational efficiency the domain only 
extends one chord length upstream of the model. (Previous investigations [16] indicate that this is sufficient 
to avoid undue upstream influence.) Since the only available experimental data for the upstream boundary 
layer is at x/c = —2.14, a turbulent boundary layer solver is used to march a solution from the leading edge, 
x/c = —4.61, and obtain a boundary layer profile at x/c = —1. In order to validate this profile, the profile 
at x/c = —2.14 was also computed; it compares well to the experimental data, see figure 3. 

In order to attenuate the outflow disturbances an exponentially stretched rectilinear grid with 100 points 
in downstream extent is patched on at the outflow and extends to x/c = 11. The increased grid spacing not 
only increases the artificial viscosity, but, more importantly, causes the FSM contribution to approach the 
RANS limit, providing a natural damping region that tends to a steady turbulent boundary layer profile. 

Both 2-D and 3-D calculations are performed. As we expected based on earlier work, [6], for the 
unforced case in which very large amplitude 2-D structures are present, secondary 3-D instabilities become 
important to the flow, as indicated by the fact that good agreement can be obtained only for 3-D calculations. 
For the steady suction case the amplitude of the 2-D structures is much reduced, as is the role of the 3-D 
structures, and good agreement is obtained using 2-D calculations. 

Except where otherwise noted, the spanwise domain was 0.16c. Periodicity was enforced on the span- 
wise end planes. Side wall effects were not accounted for in the primary simulations, although further 
investigation is underway to examine the influence of the side walls. 

The full grid with the ramping region uses 90 1x181x17 grid points in the downstream, wall normal, 
and spanwise directions, respectively. For reference, this is roughly half the downstream and wall normal 
resolution and one-tenth the spanwise resolution used in the coarse DNS [17]. 

The forcing slot was modeled by a transpiration surface with the correct mass flow rate. For unsteady 
forcing it is important to angle the jet appropriately to assure that the disturbances are introduced into the 
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Figure 3: Comparison of the experimental data and turbulent boundary layer solution at x/c = —2.14. 




/ 


Figure 4: Schematic of slot with blowing (left) and suction (right). 


shear layer, rather than blowing through the boundary layer. In the later case we have shown that the effect 
of the forcing actually decreases for high amplitude forcing [6]. For the suction case, however, wall normal 
suction was used, since it is not expected that the flow will turn to as steep an angle as the jet nozzle angle 
(see figure ). 

The initial condition for all cases was computed using CFL3D in time steady RANS mode. 
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Introduction 

Separated boundary layer flows with subsequent downstream re-attachment are very common in engineering 
applications. Over the past decade there has been extensive research on control of flow separation using 
microjet actuators. Accurate numerical solution of such flows continues to be a challenging problem mainly 
because of their non-equilibrium nature. The present work involves the use of Reynolds Averaged Navier 
Stokes (RANS) and Detached Eddy Simulation (DES) tools in the numerical validation of the ‘Case 3’ 
problem, which involves high Reynolds number flow over a two dimensional bump. In this paper, the 
RANS results will be presented and the DES is expected to be completed by the time of the workshop. 


Solution Methodology 

For the numerical simulations, the compressible RANS solver OVERFLOW [1] is used. The one equation 
Spalart-Allmaras (S-A) [2] turbulence model is used for the computations. The inviscid fluxes are computed 
using a third order upwind scheme with the Roe riemann solver. Alternatively, second or higher order 
central differencing can also be used. The viscous terms are computed using second order differences. Time 
accurate computations are performed using the implicit second order backwards difference scheme with 
Newton sub-iterations. For simulations with turbulence model modifications and detached eddy simulations, 
the TURNS code [3] will be used. The TURNS code solves the compressible RANS equations with the S-A 
turbulence model along with third order upwinding for the inviscid terms and second order differencing for 
the viscous terms. The results from the TURNS code are seen to be identical compared to the OVERFLOW 
code when both programs are run with the same input conditions. 




Figure 1: View of Grid 2 (Zoomed up view of slot region on the right). 
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Figure 2: Wall normal spacing. 



Figure 3: Upstream Velocity (x/c=-2.14) 


The 4-Block structured grids (Fig. 1 shows the coarse grid) provided by the workshop organizers are 
used in the simulations. The fine grid (Grid 1) has 210060 points and the coarse grid (Grid 2) uses every 
other point of the fine grid. Fig. 2 shows that the normal spacing at the wall is less than 1 wall unit for 
both grids. Also, Grid 1 has around 150 normal points in the upstream boundary layer - hence these grids 
appear to be sufficient for near- wall RANS calculations. Alternate hyperbolic grids with varying grid point 
distributions in the shear layer were also used, but will not be reported since the solutions are fairly grid 
independent when the wall normal spacing is less than 2.5 wall units. 

The effect of the steady/oscillatory suction/blowing is modeled as a surface boundary condition on the 
bottom of the bell-shaped chamber (zone 2). The boundary condition procedure is as follows: The density 
and pressure are extrapolated from the interior of the flow field and the momentum components are specified 
in accordance with the known mass-flow rate. Alternately, extrapolating the stagnation enthalpy instead of 
density did not make a significant difference. For the no-flow case, this boundary is treated like an inviscid 
wall. 

The upstream boundary (x/c=-6.39) was treated as a subsonic inflow condition with characteristic ex- 
trapolation. The resulting inlet profile at x/c=-2.14 is shown in fig. 3. The downstream (x/c=4.0) boundary 
conditions are set by specifying the back pressure and extrapolating all the other primitive variables. 


Implementation and Case Specific Details 

The flow in question is a non-equilibrium flow and hence separation and re-attachment could involve un- 
steady vortex shedding. Reynolds averaging can be expected to cause errors in the final solution because of 
the difficulty in characterizing the different scales that contribute to mixing. However, before using a higher 
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Figure 4: Comparison of Turbulence model and Inviscid differencing with 2D computations using Grid 2 





Figure 5: Reynolds stress prediction with the SA models and SST model at x/c locations: -2.14,1.0 and 1.2 
(from left to right) with Grid 2 


end model like hybrid RANS-LES, it is necessary to ensure that the RANS solution is not affected by other 
procedures within the framework. 

The one equation Spalart-Allmaras and two equation Menter SST [4] models are used. Fig. 4 shows the 
results of a 2D simulation. The difference between upwinding and central differencing is found to be very 
small. Both turbulence models under-predict the suction peak and the separation region and extent is poorly 
predicted. The effect of grid adaptation and refinement in the separated region and around the shear layer 
did not result in any improved predictions. 

The S-A model [2] uses a production term given by: 

- - v 

P = c b xSv u with S = S + (1) 

Three different choices of S (from the literature) were used: Sb ase u ne = \w\,S strain = \D\ and S s t ra in,vor — 
\uj\ + 2ram(0, \D\ — |u;|), where \uj\ and \D\ are the magnitudes of the vorticity and strain tensors respec- 
tively. Fig. 5 shows varied predictions of the Reynolds stress for the no-flow case. The baseline version is 
used in all the computations that follow. 

The experiment is made nominally two dimensional by the introduction of end-plates and hence a full 3D 
simulation should include them. However, solely for comparison purposes, a 3D simulation was performed 
with a viscous wall on one side and assuming symmetry at the center-span. Grid 2 was extruded in the span- 
wise direction consisting of 51 stations with a spacing of y + — 2 at the viscous wall. The 3D simulations 
seem to predict the suction peak more accurately. 

For the steady suction case, a mass flux (pv) / (poo^oo) = 0.01235 is specified at the bottom of the 
bell-shaped chamber. For the oscillatory case, a mass flux (p(t)v(t))/(p X) U 00 ) = 0.01235 sin(27 rft) is 
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Figure 6: Comparison of 2D and 3D runs for No-flow (left) and Suction (right) cases 



Figure 7: Computed Mass flow ratio across the slot on the hump surface (612 steps per cycle for oscillatory 
case). 


specified. The computed mass flow ratio (m)/(r/io 00 f/ 00 74rea s / 0 t) at the slot is shown in fig. 7. For the 
the oscillatory case, time accurate computations were performed with 3 different time-steps using 306, 612 
and 1224 time-steps per cycle with 5 Newton sub-iterations. The time-average for the latter two cases is 
seen to be identical. The baseline and controlled pressure and skin friction coefficients (time averaged for 
oscillatory case) are shown in fig. 8 and corresponding stream-line plots are shown in fig. 9. 

Overall, it is seen that discrepancies exist between the RANS results and experiments for both the 
controlled and uncontrolled cases. These could result from a number of factors: 

- Inadequacy of the RANS closure (including the Boussinesq assumption). 

- 3D and blockage effects. 

- Boundary condition specification (especially for the control cases). 

- To an extent, measurement conditions and uncertainties. 
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Figure 8: Pressure and skin friction co-efficients using S-A model and Roe Upwinding, Grid 2 
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Figure 9: Streamlines for baseline (left) and steady suction cases (right), Grid 2 
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CASE 3: RANS AND URANS APPLICATION WITH CFL3D 


C. L. Rumsey 


Computational Modeling & Simulation Branch, NASA Langley Research Center, Hampton, VA 23681-2199 


Solution Methodology 

This case was run using CFL3D, a multi-zone Reynolds-averaged Navier-Stokes code developed at NASA 
Langley [1]. It solves the thin-layer form of the Navier-Stokes equations in each of the (selected) coordinate 
directions. It can use 1-to-l, patched, or overset grids, and employs local time step scaling, grid sequencing, 
and multigrid to accelerate convergence to steady state. In time-accurate mode, CFL3D has the option to 
employ dual-time stepping with subiterations and multigrid, and it achieves second order temporal accuracy. 

CFL3D is a finite volume method. It uses third-order upwind-biased spatial differencing on the con- 
vective and pressure terms, and second-order differencing on the viscous terms; it is globally second-order 
spatially accurate. The flux difference- splitting (FDS) method of Roe is employed to obtain fluxes at the 
cell faces. It is advanced in time with an implicit three-factor approximate factorization method. 


Model Description 

For this test case, three different turbulence models were used. The first is the one-equation Spalart-Allmaras 
model (SA) [2], the second is the two-equation shear-stress transport model of Menter (SST) [3, 4], and the 
third is an explicit algebraic stress model (EASM-ko) in k-u form [5]. The first two models are both linear 
eddy-viscosity models that make use of the Boussinesq eddy-viscosity hypothesis, whereas the EASM-ko is 
a nonlinear model. The equations describing these three models can be found in their respective references. 

In CFL3D, the models are implemented uncoupled from the mean-flow equations. They are solved using 
a three-factor implicit approximate factorization approach. The advection terms are discretized with first- 
order upwind differencing. The production source term is solved explicitly, while the advection, destruction, 
and diffusion terms are treated implicitly. For EASM-ko, the nonlinear terms are added to the Navier-Stokes 
equations explicitly. 


Implementation and Case Specific Details 

Three flow conditions are computed over the hump model: (1) no flow control, (2) steady suction flow 
control, and (3) oscillatory synthetic jet flow control. The control is applied near x/c — 0.65 on the back 
side of the hump, near where the flow separates in the un-controlled state. The freestream Mach number 
is M = 0.1, and the Reynolds number is approximately Re = 936,000 per hump chordlength. For the 
oscillatory case, the oscillation frequency is 138.5 Hz. 

All computations performed for this case were 2-D. The grid used was the supplied 2-D structured grid 
number 1 (which contains 4-zones connected in a 1-to-l fashion, and approximately 210,000 grid points), 
as well as a medium-level grid made from the fine grid by extracting every-other point in each coordinate 
direction (2-D structured grid number 2). The SA model was solved on both the fine and medium grids, 
whereas the SST and EASM-ko models were only solved on the medium grid. For the oscillatory case, only 
the S A model on the fine grid was used. 

For the no-flow-control and steady suction cases, CFL3D was run in steady-state mode, utilizing local 
time-stepping to accelerate convergence. For the oscillatory case, the time step chosen yielded 360 time 
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Figure 1 : Convergence of subiteration residual during time-accurate oscillatory computation, S A model, 
fine grid. 


steps per cycle of the forcing frequency, and 5 subiterations were employed per time step. For this case, this 
number of subiterations was enough to reduce the L 2 -norm of the subiteration density residual by about 2 
orders of magnitude. See Fig. 1, which shows subiteration residual over the course of 7 time steps (with 5 
subiterations per time step) during part of the unsteady cycle. 

The boundary conditions were as follows. At the floor and hump surfaces, as well as at the side walls 
inside the cavity, solid wall adiabatic boundary conditions were applied. At the front of the grid, which 
extended to x/c — -6.39, far-held Riemann-type boundary conditions were applied. At the downstream 
boundary (at x/c — 4.0) the pressure was set at p/p re f = 0.99947, and all other quantities were extrapolated 
from the interior of the domain. This back pressure was determined via trial and error, to achieve approxi- 
mately the correct inflow conditions. At the bottom of the cavity, different boundary conditions were applied 
depending on the case. For the no-how-control case, this wall was treated as an inviscid wall. For the steady 
suction case, the boundary condition set the velocity components as follows: 

[7 = 0 V= (pV) set /p (1) 

and (pV) set was chosen in order to achieve the equivalent of a mass how of 0.01518 kg/s through a 23-inch 
span. The value used turned out to be (pV) set = -.001248p re fa re f, where a re f is the reference speed of 
sound. For the oscillatory how case, the boundary condition set the velocity components as follows: 

[7 = 0 V= [(pV) ma Jp\cos(27rFt) (2) 

where F is the frequency and t is the time, and ( pV) max was chosen in order to achieve a maximum velocity 
magnitude near to the target of approximately 26.6 m/s out of the slot during the cycle. Fig. 2 shows the 
velocity magnitude at a position near the center of the slot exit as a function of nondimensional time, on a 
line even with the hump surface, over the course of one cycle of oscillation. The value of ( pV) max used to 
achieve this condition was ( pV) max = 0.001p re fa re f. Along with the above boundary conditions given by 
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Figure 2: Velocity magnitude at the slot exit for oscillatory case, SA model, fine grid. 


Eqs. 1 and 2, the density and pressure at the bottom of the cavity are extrapolated from the interior of the 
domain. 

The top tunnel wall was treated as an inviscid wall for all of the computations submitted to the workshop. 
However, the effect of making the top wall a viscous wall was also investigated (using a different grid with 
appropriate finer normal spacing near the top wall). Resulting surface pressure coefficients are shown in 
Fig. 3. Using viscous top wall lowers the peak C p levels over the center of the hump, in better agreement 
with experiment. However, the effect does not fully account for the difference between CFD and experiment. 
Also, the viscous upper wall does not impact the C p levels in the separated region to any significant degree. 

As a nonlinear model, EASM-ko can do a better job predicting the turbulent normal stresses than linear 
models. This can be seen in Fig. 4, which shows results predicted by the three models at x/c — —2.14. 
The linear models show no perceptible difference between u'u' and v'v', whereas EASM-ko does predict a 
normal stress difference. The u f u f from EASM-ko is in better agreement with the experiment, although its 
peak near the wall is still too low. There were no measurements of v f v f . 

Examples of the effects of grid and turbulence model can be seen in Fig. 5. This figure shows streamwise 
velocity profiles in the separated region at x/c — 0.8 for the suction case, using SA on both the fine and 
medium grids, and using SST and EASM-ko on the medium grid. There is essentially no difference between 
the SA results on the fine and medium grids, and the two-equation models both predict a qualitatively 
different profile than SA at this location. 
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Figure 5: Velocity profiles at x/c = 0.8 for the suction case. 
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CASE 3: TWO-DIMENSIONAL SIMULATION OF FLOW OVER A 
HUMP MODEL WITH MENTER’S SST MODEL 


V. Katam, H. Chen, L. Huang, V. Parimi, R. P. LeBeau, and P .G. Huang 

Department of Mechanical Engineering, University of Kentucky, Lexington, KY 40506-0108 

Introduction 

Flow control of separation is an ongoing topic of investigation. This investigation of the 
phenomenon is based on our CFD experience working on flow control approaches such as 
morphing wings[l] and blowing and suction jets [2,3]. The latter studies were performed with the 
CFD code GHOST. It is this code and its techniques that we have applied to solve the current 
flow control problem. 

Solution Methodology 

All present computations were performed with the CFD code, GHOST. GHOST is an in-house 
CFD code developed at University of Kentucky by P. G. Huang. The code is based on a finite 
volume structured formulation with chimera overset grids. The QUICK and TVD schemes are 
applied to discretize the convective terms in the momentum and turbulence equations, 
respectively; the central difference scheme is used for the diffusive terms and the second order 
upwind time discretization is employed for the temporal terms. This code has been tested 
extensively and is routinely used for turbulence model validation [4,5,6]. The turbulence model 
used in the present computation is Menter's SST two equation model [7], which provides 
excellent predictive capability for flows with separation [8]. The multi -block and chimera 
features of the code allow the use of fine gird patches near the jet entrance and in regions of 
highly active flow. The code also employs MPI parallelization to allow different computational 
zones to be solved on different processors. The computations were performed on multiple 
commodity PC clusters housed with our research group at the University of Kentucky. 

Model Description 

The turbulence model used is Menter’s SST two equation model[7]. 

Implementation and case specific details 

Grid #1 (Dense grid) and Grid #2 (Coarse grid) provided on the workshop website were used for 
the computation. 

Specifying the boundary conditions proved to be a key question in this case. The 
upstream portion of the grid until x/c = -2.14 was removed so that the measured u velocity profile 
and streamwise turbulent intensity values were given as the inlet boundary condition. Several 
boundary conditions such as inviscid wall, freestream, and no-slip (viscous) wall were tested for 
the upper and the bottom boundary of the slot cavity. The best results were obtained for the 
boundary condition of a viscous wall on top and a viscous wall at the bottom of the suction cavity 
for the no flow case (Figure [1]). For steady suction cases, the boundaries remained the same 
except the bottom cavity boundary had a fixed velocity producing the workshop-provided mass 
flow rate. 
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A simulation with steady suction and a simulation with no suction were performed on 
each grid. In each case, a steady-state and a time-dependent simulation were calculated. The 
steady-state runs were done for a sufficient number of iterations until the flow data has converged 
to a constant solution. Turbulence statistics were obtained from the time-dependent simulation; 
all other data was taken from the steady-state results. Since the time-dependent simulations were 
performed only to obtain turbulence data, we used only a single timestep value for each case, 
corresponding to 0.001 seconds. 



X 


Figure 1: Boundary conditions used for the computation. All walls are viscous. 
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CASE 3: Two-Dimensional Flow Control Analysis on the Hump Model 

Sally A. Viken 

Flow Physics and Control Branch, NASA Langley Research Center, Hampton, VA 23681-2199 

Introduction 

Computational analyses have been conducted on the Wall-mounted Glauert-Goldschmied 
type body (“hump” model) with the Full Unstructured Navier-Stokes 2-D (FUN2D) flow solver 
developed at NASA LaRC. This investigation uses the time-accurate Reynolds-averaged Navier- 
Stokes (RANS) approach to predict aerodynamic performance of the active flow control 
experimental database for the hump model. The workshop is designed to assess the current 
capabilities of different classes of turbulent flow solution methodologies, such as RANS, to 
predict flow fields induced by synthetic jets and separation control geometries. The hump model 
being studied is geometrically similar to that previously tested both experimentally and 
computationally at NASA LaRC [ref. 1 and 2, respectively]. 

Solution Methodology 

The FUN2D flow solver is a node based, implicit, upwind flow solver used for 
computing flows around airfoil configurations discretized with an unstructured grid [ref. 3]. The 
governing equations (provided below) are the time-dependent RANS equations in conservation- 
law form, which are integrated in time to obtain a steady state solution. The inviscid fluxes are 
obtained on the faces of each control volume by using the flux-difference-splitting (FDS) 
technique of Roe [ref. 4]. A node-based algorithm is used in which the variables are stored at the 
vertices of the mesh and the equations are solved on non-overlapping control volumes 
surrounding each node. The viscous terms are evaluated with a finite-volume formulation that 
results in a central-difference-type scheme. The Spalart-Allmaras (SA) turbulence model is used 
in this investigation and all computations assume fully turbulent flow [ref. 5]. 

A two level iteration is used to achieve convergence of the discrete algebraic equations at 
each time-step. The outer iteration is a modified Newton method and employs a first-order Van 
Leer Jacobian (LHS) [ref. 6] driving the second-order residual vector (RHS). The inner iteration 
employs a red-black Gauss Siedel point-implicit algorithm to solve the equations at each step of 
the outer iteration. Twenty inner sub-iterations were used in all cases in this study. The SA 
turbulence model equations are weakly coupled to the hydrodynamic equation via the outer loop 
iteration. 

For steady state computations, the solution is driven to convergence using an Euler 
implicit advancement in pseudo-time. For time-dependent computations, the solution is 
discretized in physical time with the second-order backwards differentiation formulae (BDF), 
while pseudo-time iterations are again employed to relax the equations [ref. 2, 7]. The linear 
system of equations resulting from either formulation is iteratively solved with a point-implicit 
procedure. 

Governing Equations 

The governing equations, the Reynolds-averaged Navier-Stokes equations, are written in 
conservative-law form to relate the rate of change of mass, momentum, and energy in a control 
volume of area A to the fluxes of these quantities through surface of the control volume. The 
non-dimensional equations presented in vector form are as follows: 

A — + <J)Fi-n dl - <j)fVnd/=0 

Xi xi 
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where n is the outward-pointing unit normal to the surface of the control volume 3Q. The 
solution vector Q containing the dependent variables is defined by 



The inviscid and viscous flux vectors through the surface of the control volume 3Q, defined as F i 
and F v? respectively, are given by 
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The normal and shear stress terms and heat conduction terms are defined by 
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The equation of state for a perfect gas is used to define the pressure p 
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and the laminar viscosity u is determined through Sutherland’s law 
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where C* = : — is Sutherland’s constant divided by a free stream reference temperature which 

536.4 

is defined to be 536.4° Rankine for Case Study #3. 

The eddy viscosity, |i„ is obtained by the turbulence closure model developed by Spalart 
and Allmaras [ref 5], Again, the turbulent viscosity equation is solved separately from the flow 
equations at each time step, using the Euler implicit time-stepping scheme, resulting in a loosely 
coupled solution process. 

Implementation and Case Specific Details 

Unstructured Grids 

Two-dimensional unstructured grids were generated for Case Study #3 analyses and 
provided at the website. The coordinates used for the unstructured grids were the 2-D non- 
dimensionalized theoretical model coordinates along with the tunnel geometry. These grids were 
generated with advancing front type point placement with iterative local re-meshing for grid 
quality improvement [ref. 8,9]. For the internal flow analyses conducted with the FUN2D code 
that will be presented at the workshop, the website unstructured Grid #1 and #2 were used. The 
forward extent of these grids is longer than the actual splitter plate length used in the wind tunnel 
experiment. These grids ran from -6.39c ahead of the leading edge of the model to 4.0c behind 
the leading edge of the model. This grid extent of -6.39c was chosen because it was found in 
preliminary CFD tests to yield a "run" long enough so that the computed boundary layer 
thickness approximately matched that of experimental data at x/c = -2.14 [figure 1]. The grid 
height is y/c = 0.90905, which corresponds with the actual height of 15.032 inches from the 
splitter plate to the ceiling of the tunnel test section. In addition, these grids were generated with 
the internal cavity modeled. The fine grid (website unstructured Grid #1) has 123703 nodes, 
368476 faces, and 247404 cells. The minimum spacing at the viscous walls was set to be 
approximately 8.e-06. This minimum spacing yields a y+ value less than 1. The coarse grid 
(website unstructured Grid #2) has 57152 nodes, 169689 faces, and 114302 cells. The minimum 
spacing at the viscous walls was set to be approximately 1.6e-05. The upper grid boundary (test 
section ceiling) was set up to use inviscid-type grid spacing. 

Initial calculations were conducted on the hump model using the grids described above for 
the baseline no-flow control condition. The experimental data for the baseline case showed a 
discontinuity in the Cp distribution at x/c = 0.48, which was not observed in the computational 
results. An unstructured grid was generated using the Quality Assurance (QA) coordinates 
provided at the website followed by further computational analyses. This grid has 122790 nodes, 
365737 faces, and 245578 cells, with the minimum wall spacing yielding a y+ value less than 1. 
Computational results for the baseline no-flow condition analyzed with the QA coordinate grid 
and theoretical coordinate grid are compared with the experimental data [figure 2], Computations 
with the QA grid show the sensitivity of the solution to the discontinuity in the hump model 
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surface geometry similar to that found experimentally. Outside of this localized discontinuity in 
the Cp distribution, the CFD solutions obtained with the QA coordinate grid and the theoretical 
coordinate grid compared well together. Grid studies were also conducted with the upper grid 
boundary defined to have viscous wall spacing, however only a slight increase in the flow 
acceleration over the hump model was achieved. All computations that will be presented at the 
workshop were run using the grids generated with the theoretical coordinates (unstructured Grids 
#1 and #2) that have inviscid-type grid spacing on the upper grid boundary. 

Boundary Conditions 

The boundary conditions on the hump model, the internal cavity walls, and the tunnel 
floor corresponded to no-slip between the fluid and the solid boundary at their interface, with a 
constant temperature wall that was set to the adiabatic wall temperature T aw . The tunnel ceiling 
was treated as an inviscid surface. The internal "actuator" boundary at the bottom of the hump 
cavity was also treated as an inviscid surface for the baseline cases with no control. For the 
control cases, the "actuator" boundary condition corresponded to pV(q, t| = 0, t) = 
Amplitude *f(^)*cos (cot) where f(q) = 1 (tophat distribution). The mass flux through the 
“actuator” boundary of the internal cavity was adjusted to obtain a peak velocity near 26.6 m/sec 
at the slot exit during the blowing part of the cycle [(pV) max = 0.001(p re fa re f) was used to obtain 
this condition]. For inflow conditions, temperature was also specified on the "actuator" boundary 
from the experimental data. To obtain the boundary conditions at the tunnel inlet, the flow was 
assumed to be both inviscid and isentropic in this region so that quantities for the computation of 
the flux along the inflow boundary were obtained from two locally 1-D Riemann invariants. The 
Riemann invariants were considered constant along characteristics defined normal to the inflow 
boundary. At the downstream boundary, a back pressure of 0.99947 times reference pressure was 
specified in order to approximate the upstream conditions at the tunnel inlet. 

Time-Step 

The baseline and steady suction cases were run non-time-accurate, and achieved steady state 
convergence. The oscillatory suction/blowing case was run time-accurately using the second- 
order accurate Backward Differentiation Formulae (BDF) scheme. The assumption was made that 
the excitation frequency introduced at the “actuator” boundary was a perfect sine wave. The 
computational oscillatory cases were run time-accurately at least 20 shedding cycles to set up the 
flow field and then the Cp’s were averaged over 20 cycles. 

Test Conditions 

Computational flow analyses were conducted on the hump model for M = 0.1, at Re = 
936000. The two required test conditions for Case Study #3 were conducted along with the 
optional condition; the no-flow through the span slot, the suction case with suction rate of 

0. 01518. kg/sec through the slot, and the zero-net-mass-flux oscillatory suction/blowing case 
(frequency = 138.5 Hz, and peak velocity out of slot during blowing part of cycle = 26.6 m/s). 
Only two-dimensional computations were performed for these studies. 

An example of the grid density effects are shown in figure 3, where the streamwise 
velocity profiles in the separated flow region at x/c = 0.8 for the suction case can be observed. 
The profiles for the fine and coarse level grids show that the grids are sufficiently resolved for the 
case analyzed. 
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Figure 1. Velocity profile at x/c = -2.14 for 
baseline case (1VL = 0.1; Re = 936000). 



Figure 2. Baseline computational results 
versus experimental data (1VL = 0.1; Re = 
936000). 



Figure 3. Velocity profiles at x/c = 0.8 for 
the suction case (1VL = 0.1; Re = 936000). 
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Case 3: Separated Flows — An Assessment of the Predictive Capability of Five 

Common Turbulence Models 


R.E. Spall, W.F. Phillips, and N.R. Alley 
Department of Mechanical and Aerospace Engineering 
Utah State University, Logan, UT 84322-4130 


Introduction 

The results presented herein were generated using the commercial CFD solver Fluent, which 
makes available to the user a relatively large number of turbulence model options. Fluent is 
widely used in both industry and academia, and it was of interest to the authors to assess the appli- 
cability of the available models in predicting separated flows. Our interest is also derived, in part, 
from our use of the code in the design of high lift, unmanned aerial vehicles (UAV’s) in which the 
ability to predict flow separation is crucial. 

Solution Methodology 

The results were computed using the commercial CFD solver Fluent (version 6.18). Fluent 
employs a pressure-based finite volume solution procedure to solve the governing equations on 
unstructured grids. Pressure-velocity coupling was accomplished using the SIMPLEC procedure. 
Second-order central differencing was used for the viscous terms in the transport equations. Inter- 
polation to cell faces for the convection terms was performed using a third-order QUICK scheme. 
In all cases, a steady, two-dimensional formulation was employed. Solutions obtained using a seg- 
regated solver were considered converged when residuals for each of the equations (based on an 
L2 norm) were reduced by a minimum of five orders of magnitude. Additional iterations were 
then performed to confirm iterative convergence. 

Model Description 

Five different turbulence models available in the commercial CFD solver Fluent were used to 
solve test case 3. These models included 1) k- £ , 2) k-ti), 3) Reynolds stress transport, 4) Spal- 
art-Almaras, and 5) shear stress transport (SST). In all cases, the near wall mesh within the chan- 
nel region along the lower wall was sufficiently fine so that wall functions were never 
implemented. A brief description of each model appears below, with an emphasis on describing 
methodology that often varies within each class of model. 

1) The k - £ model 

The standard k-e model of Launder and Spalding [1] was employed, with low-Reynolds num- 
ber modifications for the near-wall modeling, which combine a two-layer model with wall func- 
tions. In the present work, the near wall mesh was fine enough to resolve the viscous sublayer; 
consequently the two-layer model was used and wall functions were never implemented. In par- 
ticular, if Re y < Re* (where Re* = 200 ) the one-equation model of Wolfstein [2] is employed. 
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In this region, the equation for k is retained, but the turbulent viscosity is obtained from 

\i t • = p C^l^Jk where the length scale is given as = yc { ( 1 - e Rey/A,i ) [3], The turbulent vis- 
cosity determined above is then smoothly blended with the high Reynolds number (i f 0 obtained 

in the outer region. The dissipation rate in the near wall region is also specified algebraically as 

3/2 

£ = k /l £ where the length scale, l e , is computed using the same relation used in the specifi- 
cation of the turbulent viscosity, although with a different value of the constant A . A blending 
function is used to ensure a smooth transition between £ specified algebraically in the inner 
region, and £ computed via the transport equation in the outer region. The blending is of the form 

v t = 4^.0 + a (1) 

where. 


X 


1 

2 


1 + tanh 


(Re y -Re* y ) 

A 


( 2 ) 


2) The k - to model 

A low-Reynolds number version of the Wilcox k - (0 model [4] was employed. In particular, a 
low Reynolds number correction to the turbulent viscosity, |l f = a*p£/co, is determined as: 


a = a 


Vp + Re/Re^ 
1 +Re/Re kJ 


( 3 ) 


where Re t = (pfc)/(|i(0) , R k = 6, 0Cq = (3 t -/3 , = 1 , and (3- = 0.072. A low-Reynolds 

number correction is also applied to the production of 0) term, G m = OUi)G k /k, as: 


a = — 

* 

a 


a~Ya 0 + Re/R m 


1 +Re/R a 


( 4 ) 


where = 0.52 and R a = 2.95 . In terms of boundary conditions, the asymptotic value of (0 
at the wall is specified as: 

to w = p(«*) 2 (0 + /p. (5) 

The value of (0 + is determined from: 

6 


0) = mm 


2500, 


Pl(/) 2 - 


( 6 ) 


where f/ = 0.09. 


3) Reynolds stress transport model 

The RST model also uses a two-layer formulation for the near wall, similar to that used for the 
k - £ model. In this case, modifications to the linear pressure-strain model follow the methodol- 
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°gy given in Launder and Shima [5]. These modifications involve specification of the constants 
multiplying the slow and rapid pressure strain terms, and the wall-reflection term, in terms of 
Reynolds stress invariants and the turbulent Reynolds number. Diffusive transport is modeled 
using a scalar turbulent diffusivity [6], given as: 


D 


T.ij 


_d 

dx 


k\°k dx k J 


(7) 


4) Spalart-Allmaras 

The Spalart-Allmaras model [7] solves a transport equation for a variable that is a modified form 
of the turbulent kinematic viscosity. For the results presented, the deformation tensor, 5, appear- 
ing in the production term follows the original model proposed by Spalart and Allmaras, which is 

based on the magnitude of the vorticity S = , where kl- is the mean rate-of-rotation 

tensor. Fluent also incorporates a modified definition of S which includes measures of both rota- 
tion and strain tensors in its definition; however, this modification was not used. 


5) Shear Stress Transport 

Fluent provides a low-Reynolds number version of the SST model, which was employed for the 
results presented. In particular, the turbulent viscosity was determined as: 

1 


ok 

fir = “ 


(0 


max 


I QF 2 


a * ctjCO 


( 8 ) 


where a* is as defined for the k- 0) model and Q = . In addition, the turbulent 

Prandtl numbers are defined as: 
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= 


<*(0 = 


^ G k, 1 + (1 _ 2 

1 

7V°(o,l + (!- F i ) /( V 2 


where F l and F 2 are blending functions. In addition, <3 k j anva kj ^ 2 
turbulent Prandtl numbers, respectively. (Similarly for G 0) , and o 0) 2 .) 


(9) 

( 10 ) 

and G. t are constant inner and outer 


Implementation 

The inlet n-velocity boundary condition was specified using the experimental profile available on 
the web site. The v- velocity was set to zero. Turbulence quantities were derived from the approxi- 
mate inlet turbulence intensity (7) of 0.09% (given on web site) and a turbulence length scale (/) 
which was defined as 0.4 8 0 99 (where 8 0 99 is the boundary layer thickness at the inlet). In partic- 

2 

ular, the inlet turbulence kinetic energy was then set to k = 1.5 (ul) where u is the mean flow 
velocity. The dissipation rate was computed as £ = C^' 75 k L5 // . For the Spalart-Allmaras model. 
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the inlet modified viscosity was specified as v 


3- 

-ull. When implementing the k- 0) model, 


1/2 1/4 

the inlet specific dissipation rate was computed from G) = k /( l) . Finally, for the RST 

model the inlet Reynolds stresses were specified as u-uj = 0 and = (2/3 )k, where k 

was computed as described above. In all cases, zero normal derivative boundary conditions were 
used at the outflow plane. No-slip conditions were implemented on both the upper and lower 
walls. The model-specific approaches to implementing wall boundary conditions were described 
in the section on Model Description. (Slip wall conditions were also used on the upper wall, with 
little change in the resulting solutions. However, the results presented are for the no-slip condi- 
tions.) 

Solutions were computed on two different grids, neither of which were those provided on 
the web site. The coarser grid consisted of 38,720 quadrilateral cells within the channel; the finer 
grid consisted of 85,760 cells. Grid points were clustered toward the wall so that in all cases 

y + < 0.5 for the lower wall adjacent cells. This was achieved by specifying the distance of the 
first grid point above the lower wall as y = 0.000018c . The grid spacing along the upper wall 

was defined such that y + ~ 50 , so that wall functions were used. For the suction case, the suction 
chamber was included in the model and contained 2,169 tetrahedral cells. No grid refinement was 
performed within the suction chamber. The grids (either fine or coarse) within the channel section 
were identical for both the no-suction and suction cases. A section of the finer grid in the region of 
the suction slot is shown in Fig. 1. 
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Figure 1 . Mesh for Case 3 with flow control in the region of the suction slot. The 
grid consists of 85,760 quadrilateral cells in the channel region, and 2,169 tetra- 
hedral cells in the suction chamber. 
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CASE 3: Active Control for Separated Flows 


A. Shmilovich, Y. Yadlin, and R.W. Clark 


Aerodynamic Technology, Flight Sciences/Advanced Design, Phantom Works, The Boeing Co., Huntington 
Beach, CA 92647 


Introduction 

A numerical procedure has been developed for the analysis of unsteady flows within the framework of 
active control. Computations have been performed for flows past airfoils, wings and engines, ranging 
from nearly static ambient conditions and up to transonic freestream velocities. The analysis tool has been 
extensively used to develop new and practical flow control approaches. Some of these applications 
include vortex control (by reducing core strength or by introduction of perturbations to destabilize the 
vortex) and separation control to enhance high lift performance over a wide range of flow conditions. 


Solution Methodology 

The numerical tool is a modified OVERFLOW code originally developed by NASA [1]. Overflow uses 
the Reynolds Averaged Navier Stokes formulation in conjunction with field turbulence models. A special 
module has been used for applying the time- varying boundary conditions. Various signal shapes can be 
used in conjunction with arbitrary stagnation properties for the general description of the jet. 


Implementation and Case Specific Details 

Several aspects of computations are described bellow. 

Grid: The point-match grid provided by the workshop organizers has been used for these simulations with 
a slight modification. A narrow region of grid overlap has been introduced by extending the grid of the 
actuator into the wind-tunnel channel grid. 

Turbulence models: SA and SST models have been used. Results in terms of pressure distributions and 
off surface streamwise velocity were very similar. 

Boundary conditions: The program uses the mass flow rate, area and the stagnation pressure and 
temperature to define the velocity across the boundary. The orientation of the flux vector can also be 
prescribed and jet pulsation is obtained from the forcing frequency. General signal shapes ranging from 
sinusoidal to step-function are defined by a set of analytical functions. A sine wave representation of the 
jet is being employed in the oscillatory case. 

Time steps: The time-accurate calculation for the pulsed jet case uses 800 steps per cycle. The calculation 
starts with a steady-state solution obtained for the flow in the absence of any excitation. Asymptotic 
quasi-state is achieved after approximately 20 cycles based on the normal force on the bump. 

With respect to the pulsating actuation case it is instructive do discuss the issue of flow separation and 
reattachment. Fig. 1 shows the off-surface streamwise component of the velocity along the bump. The 
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instantaneous distributions at equal time intervals over one activation period (0.00722 seconds) show the 
transitory nature of the separation bubble. The onset of flow separation is directly affected by the 
intermittent ejection and suction at the orifice and the entrainment of the surrounding fluid. The size of 
the separation region and the movement of the reattachment point are being influenced by the oncoming 
flow and the impact occurs on a larger scale. 


0 . 08 ,— 



Figure 1 : Extent of flow separation on the bump based on instantaneous off-surface 
streamwise component of the velocity (pulsed excitation case) 


References 

[1] Buning, P.G. et al, “OVERFLOW User’s Manual, Version 1.8j,” NASA Langley Research Center, 
Hampton V A, 1999. 


3.12.2 




CASE 3: Investigation of a boundary condition model 
for the simulation of controlled flows 


C. Marongiu 1 , G.Iaccarino 2 , P. Catalano 1 , M. Amato 1 

r CIRA Italian Center for Aerospace Research Via Maiorise, Capua (CE) 81043 - Italy 
2 CTR Center for Turbulence Research Stanford University, Stanford CA 94305 - USA 

Introduction 

The aim of this work is to investigate the use of a boundary condition for the simulation of flows 
controlled by steady-unsteady mass injection/ suction devices. These actuators usually consist of a 
deforming membrane oscillating inside a cavity. The simulation of the whole cavity would require a 
very complicated and time-consuming moving-mesh calculation. Therefore the use of a boundary 
condition applied at the cavity exit, the interface between the external and the internal domain, is an 
attractive strategy. This work is a first step towards investigating the limits and the benefits of the 
approach. 

Results obtained by two unsteady RANS codes, the CIRA U-ZEN code, and FLUENT, a 
commercial CFD code, are reported for the flow around a two-dimensional bump. Several 
turbulence models are employed including one-, two- and four-equation models. The flow control is 
obtained via suction and it has been computed using two boundary conditions, corresponding to a 
top-hat distribution on the wall and to an inclined jet. The results have been compared with the 
experimental data and numerical computations including the whole cavity. 


Computational Codes 

U-ZEN 

The CIRA U-Zen code solves the compressible RANS equations around complex aeronautical 
configurations using multiblock structured grids. The numerical discretization is based on a second 
order- cell centered finite volume method with explicit (fourth order) artificial dissipation. The 
unsteady procedure is based on the dual time stepping method where a pseudo steady-state problem 
is solved at each time step. Conventional convergence accelerators, including geometrical multigrid 
and residual smoothing, in the dual integration are used. Several turbulence models are available in U- 
ZEN [1]: for the numerical simulations presented in this work the Myong and Kasagi k — s [2], and 
the SST Menter k-co [3] turbulence model are used. 

FLUENT 

FLUENT is a commercial CFD code that solves the RANS equations on hybrid unstructured grids. 
It uses a second-order upwind discretization based on the SIMPLE pressure-velocity coupling and 
the formulation can accommodate compressible flows. Dual time stepping is used to obtain time 
accurate simulations and an algebraic multigrid technique is used to accelerate convergence within 
each time step. A multitude of turbulence models and variants are available in Fluent. In this work 
the Spalart Allmaras, SA, [4] model and Durbin’s iTf four-equation model [5] are used. The i?f model 
was implemented in FLUENT using the User Subroutines [6]. 

Flow Simulations 

Basemodel 

The basemodel case has been studied as a two-dimensional flow. Several computational grids and 
turbulence models have been used. The aim is the evaluation of the accuracy of the predictions in 
terms of grid and inflow and outflow conditions. The domain extends for 2.14c upstream and 4c 
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downstream the hump. The height is 0.91c. The external boundary conditions are specified as 
follows: 


u = from experimental data 


inflow 

v = 0 


(i) 


-iU 1 + lA, 

Pn 2 Pa> 

(VA^) 


outflow 

u and v, extrapolated 
-2- = 0.99944 

P oo 

(2) 

upper wall 

slip wall 


(3) 

lower wall 

no - slip wall 


(4) 


U-ZEN computations are carried out using the available grid on the CFDVAL2004 web site, (by 
C.Rumsey, GRID1) and the k-CD and the k-£ turbulence models. The computation with the k-co 
model has been repeated with a second grid (again available on CFDVAL2004 web site, supplied by 
C. Marongiu, GRID2). These three tests on the base model configuration have been performed with 
a turbulence level intensity equal to 0.09%. A fourth test has been performed imposing a kinetic 
energy profile together the condition expressed in (1), at the inflow section, derived by the limiting 
values assumed by the velocity fluctuations at the end of the viscous sublayer in a fully- developed 
turbulent boundary layer [7] : 


W= 0.5 u' 
v f = 0.75i/ f 


where u’ is known by the experiments. Globally, U-ZEN has given similar results between the two 
grids, under the same boundary conditions, showing a strong sensibility to the inflow conditions. In 
fact, a remarkable difference of the case with the assigned ^-profile can be seen especially in the 
boundary layer. The pressure coefficient on the hump is well predicted in the expansion zone, but 
after the separation, the results are generally inaccurate with a strong dependency on the turbulent 
model. 

The FLUENT simulations have been carried out using the same grid used before (GRID1) with and 
without the cavity. The boundary conditions used are similar to what explained before, the only 
difference regards the treatment of turbulent quantities. For the SA model a constant turbulent 
viscosity ratio (equal to 10) is specified, whereas for the four equation model the same values used 
for the k-£ model have been employed. The results are generally similar to the ones obtained with U- 
ZEN in the expansion zone, with a slight overprediction of the pressure peak at the leading edge of 
the bump (this was found to be extremely sensitive to the inlet and the upper wall boundary 
conditions). The separation zone is strongly influenced by the turbulence model, but both appear to 
underpredict the pressure level measured in the experiments. In particular, the v 2 -/ model predicts a 
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lower level of the pressure but a better overall recovery in the reattachment region downstream of 
the bump. 


Steady suction case. 

For these simulations only the k-co turbulence model has been used in U-ZEN and the SA model in 
FLUENT. The suction condition has been applied at the slot exit. The mass flow rate is constant 
along every section of the slot, while the momentum coefficient depends explicitly on the velocity 
magnitude. The simplest formulation used is to consider a top-hat velocity distribution applied along 
the d/c segment (see figure 1). The momentum coefficient will be scaled as follows: 

C>*)=C>*)~ <«) 


where the superscript ( ) + refers to the dimensionless quantities. Considering the figure 1, it is clear 
that the jet does not exit along the wall normal, but is inclined with a certain angle depending on the 
lateral wall inclinations and on the interaction with the external flow. Therefore, the problem of 
definition of a tangential velocity component arises. The second formulation consists in introducing a 
realistic direction by adding a tangential component that contributes to the C and not to the mass 
flow. Assuming a vector notation, the effect of the jet angle can be considered as follows : 

u:=ui+u'j,=u:j,+u' Jn 4kt (7) 

k = tan 2 6 


c m =P + u + Jn h + 8 

C„=2p'ui(Uky <8) 

where 6 is the angle between the wall normal and the jet direction, U_ + Jn and U^ Jt are respectively the 

normal and the tangential components of U_j . Indeed, to take into account the stagnation region at 

the beginning of the orifice, (see figure 2), a value of less than the geometrical one can be 
considered. In this simulation the following values have been used: 

9 = 43.9 ° 
ct = 3.5 10' 3 

In the figure 4, 5 and 6, the u-field component of the simulations with the whole cavity (FLUENT) , 
the wall normal jet, and the inclined one (U-ZEN) , can be seen. An improvement of the results has 
been obtained with the inclined model, although more accurate models could be defined. In fact, the 
treatment of the turbulent variables at the slot exit, here considered as a solid wall, is an open matter 
that could affect the solution quality. 

Conclusions 

From the analysis of the configuration without control, a strong sensitivity of the solution to the 
inflow and outflow conditions has been found. It’s worth to note that the upper surface has been 
considered as a solid wall with a slip boundary condition, and it is not possible to exclude an 
influence of the upper boundary layer especially in the separated region (blockage effect, [8]) . In the 
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suction case, the computations simulating the whole cavity (FLUENT) and those applying the 
boundary condition at the slot exit (U-ZEN) are in a good agreement with the experimental data 
especially if the effect of the jet direction is taken into account. This is accomplished by adding a 
tangential component; the formulation of this boundary still does not include any modification of the 
turbulent variables. This would be the scope of future investigations. 



Figure 1- Slot exit detail 



Figure 2 Cavity simulation. (Fluent) Detail on the 
stagnation region. 



Figure 3 u-contours. Suction case. Experimental data 



Figure 4 u-contours. Suction case. Cavity simulation 
(FLUENT) 
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Figure 5 u-contours. Suction case. Wall normal jet. 
(U-ZEN) 


Figure 6 u-contours. Suction case. Inclined jet. 
(U-ZEN) 


References 

[1] Catalano, P., Amato, M., “An evaluation of RAN S turbulence modeling for aerodynamics applications”. 
Aerospace Science and Technology, Vol.7, No.7, October 2003, pg 493-590 

[2] Myong, H., Kasagi, N., “A new approach to the improvement of the k-s turbulence model for wall 
bounded shear flows”, JSME Intern. J. Ser.2 33 (1) 63-72. 

[3] Menter, F., R, ‘Two equation eddy viscosity turbulence models for engineering applications ”, AIAA J. 32 
(8) (1994) 1598-1605. 

[4] Spalart, P. R., Allmaras, S., “A one-equation turbulence model for aerodynamic flows”, AIAA 92-0439 
(1992). 

[5] Durbin, P. A., “Separated Flow Computations with the k-s-v 2 model”, AIAA J. 33 (1995) 659-664. 

[6] Iaccarino, G., “Prediction of a Turbulent Separated Flow Using commercial CFD Codes”. J. Fluids 
Eng., V. 123, (2001) 815-824. 

[7] Mathieu, J., Scott, J., “An introduction to turbulent flow”, Cambridge University Press 2000 

[8] Iaccarino, G., Marongiu, C., Catalano, P., Amato, M., “BANS simulation of the separated flow 
over a bump with active control ”, Annual Research Briefs 2003, Center for Turbulence Research, 
Stanford University/ NASA Ames. 


3 . 13.5 



CASE 3: SYNTHETIC JET 

3D STEADY HALF-DOMAIN COMPUTATION USING CFD++ 


S. Shariff, P. Batten 


Metacomp Technologies , 28632-B Roadside Drive, Suite 255, Agoura Hills, CA 91301 


Introduction 

This paper briefly explains the details of one of our approaches to Case 3 of the Langley Research Center 
Workshop on synthetic jets. Here we have performed a three-dimensional simulation of half of the domain. 
For the suction case, we have prescribed the outflow pressure in the synthetic jet slot. This pressure was 
then adjusted to produce the appropriate mass flow through the slot. 


Solution Methodology 

This work was performed using CFD++, a Navier-Stokes solver for compressible and incompressible flow. 
CFD++ features a second order Total Variation Diminishing (TVD) discretization based on a multi-dimen- 
sional interpolation framework, which is also utilized for the viscous terms. A pre-conditioned HLLC 
(Harten-Lax-van Leer with Contact wave) Riemann solver is used to provide proper signal propagation 
physics while preserving positivity and satisfying the entropy condition. Further details regarding the nu- 
merical methodology can be found in Chakravarthy et al [2], Peroomian et al [4, 5], and Batten et al [1]. 

The simulations were computed with Metacomp’s 2-equation cubic n-e RANS turbulence model, which 
is discussed in Palaniswamy et al [3]. 


Implementation and Case Specific Details 

The grid was generated from the supplied geometry, and consists of 2.47 million hexahedral cells. The 
portion of the grid near the airfoil is shown in Figure 1 to give an idea of the stretching used. The grid was 
tightly stretched towards the bottom of the tunnel to directly resolve the boundary layer on the airfoil, but 
less tightly stretched at other walls. The domain extent in the streamwise direction was —2.14 < x/c < 
4.94. The shorter domain was chosen because CFD++ allows the user to impose a fully-developed Musker 
boundary layer profile as a boundary condition, and the upper-wall boundary layer thickness was specified 
at x/c = —2.14. The lower boundary layer thickness for the Musker profile was estimated using the splitter 
plate length. 

The outflow condition was set to a characteristics-based inflow/outflow condition. The equations were 
solved to the wall at the airfoil surface and lower tunnel wall, and wall functions were used for all other 
walls. 

The flow in the slot was modeled, and for the suction case, a pressure-based inflow/outflow boundary 
condition was used at the bottom of the slot. This BC imposes the pressure in the case of outflow, and 
determines the state from the interior velocity in the case of inflow. Initially, the pressure was set to an 
arbitrary value somewhat less than the cavity pressure which developed in the no-suction case. This pressure 
was then adjusted to achieve the specified mass flux of 0.01518 kg/s through the slot (which corresponds to 
0.00759 kg/s through the half-slot simulated here). This mass flux was ultimately achieved by setting the 
boundary pressure to 2100 Pa lower than the freestream, or 99.225 kPa. 
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CASE 3: COMPUTATIONS OF FLOW OVER THE HUMP MODEL 
USING HIGHER ODER METHOD WITH TURBULENCE MODELING 


P.Balakumar 

Flow Physics and Control Branch, NASA Langley Research Center, Hampton, VA 23681 


Introduction 

The flow over the two-dimensional hump model is computed by solving the RANS equations 
with k-co (SST) model. 

Solution Methodology 

The governing equations, the flow equations and the turbulent equations, are solved using the 5th 
order accurate weighted essentially non-oscillatory (WENO) scheme for space discretization and 
using explicit third order total-variation-diminishing (TVD) Runge-Kutta scheme for time 
integration. The WENO and the TVD methods and the formulas are explained in [1] and the 
application of ENO method to N-S equations is given in [2], The solution method implemented 
in this computation is described in detail in [3], 

Model used 

Standard k-co (SST) model is used and the equations and the model coefficients are described in 
[4, 5, 6], 

Implementation and Details 

The computational domain extends from x/c=-10. to 4. in the streamwise direction and extends 
from the splitter plate to the upper tunnel wall in the normal direction. The leading edge of the 
splitter plate is modeled as a super ellipse with 0.25 in. half thickness and an aspect ratio of 2. 
The leading edge is located at x/c=-5.9. C-Type grid is used around the splitter plate and a 
rectangular grid is added upstream of the leading edge as shown in Fig. 1. This grid overlaps the 
C grid and 5 th order central interpolation is used to transfer the flow variables from one grid to 
the other at the boundaries. (651*151) grid size is used around the splitter plate and the hump 
and (101*51) grid size is used in the rectangular region. 
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Figure 1 : Overlapping grid near the leading edge of the splitter plate. 


Following boundary conditions are implemented at different boundaries: 


1. 


2 . 


3. 

4. 

5. 

6 . 


At the upper wall inviscid conditions are applied. 

dp _du _di E_ v _q 

dy dy dy 

At the lower wall viscous conditions are used. 
u — v — 0 , 

T =T — T 

w free stream ’ 


( 1 ) 

( 2 ) 


and p is computed from the continuity equations. 

From the leading edge of the splitter plate to the inflow boundary symmetric conditions 
are used. 

At the inflow boundary stagnation pressure, one Riemann variable and normal velocity 
v=0 are prescribed and the second Riemann variable is solved for to obtain the other flow 
quantity. 

At the outflow boundary the pressure is specified to obtain the required Mach number 
and characteristic-type boundary conditions are implemented similar to as described in 
[7] to obtain other flow variables. 

In the suction case, boundary conditions are applied on the surface of the hump across the 
suction slot. The suction slot extends from x/c=.6541 to .6584 and across the slot normal 

mass flow rate is specified to the experimental value. A suction distribution of the form 

r \ 


(Pv)n = /max ^ 


7t(x- 


,) 


( 3 ) 


end ^ start')/ 

is used. Other forms have been tried and all of them yield the same results for a fixed 
total suction rate. The other flow quantities are obtained using the characteristic type 
boundary conditions [7]. 


Following boundary conditions are implemented for the turbulent quantities at different 
boundaries. 

1. At the inflow boundary small values are prescribed for k and p T . 
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( 4 ) 



= .009. 
fL 

2. At the outflow boundary k and to are solved for from the governing equations. Higher 
order extrapolation condition is also tried and it gives the same results. 

3. At the lower viscous wall k=0 condition is used and following exact boundary condition 
is derived for (D. 

4. In the suction case, across the suction slot linear extrapolation is used to obtain the 
turbulent quantities on the surface. 


Since the variable to becomes singular near a viscous wall, in practice a large approximate value 
is prescribed at the wall. 


® W aU 


60a, 
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where d w is the distance to the first grid point from the wall. This is an approximate boundary 
condition and when it is implemented in the higher order scheme, oscillations and convergence 

problems are encountered and the following exact condition is derived for co wall . By realizing 

/ \ 


that 



type singularity for the variable to arises because of the balance between the 


dissipation term and the viscous diffusion term in the (0 equation, the singularity is removed by 
rewriting the variable to as 


C 

0)-—r0) 1 , 

y n 


( 6 ) 


where y n is the normal distance to the wall, C is a constant and CQi is the new variable which is 
now regular near the wall. When this is substituted into the to equation the following equation is 
obtained for coi, which is similar to the to equation except for the source term. 
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The value of coi at the wall becomes 


2 


= tin 

M*, 


Here the variable to is nondimensionalised by 



( 8 ) 


(9) 


Hence the procedure is to use the C0i equation for the first few points near the wall and switch to 
the to equation away from the wall. In these computations gl>i equation is solved for the first ten 
points near the wall. Figure 2 shows the distribution of k, co and (Oi near the wall and it is seen 
that this technique resolves the viscous layer smoothly even though co is infinite at the wall. 
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Figure 2: Variation of k, co and (Oi near the wall at x/c=0.4. 

Figure 3 shows the contours of the U velocity near the leading edge region and for the entire 
computational domain. Near the leading edge region the flow separates and forms a small 
separation bubble. In the no flow case, the solution converged for the flow equations and the 
turbulent equations very well. In the suction case, the maximum residual in the co equation near 
the suction slot edges converged only by three orders. This may be due to the boundary 
conditions used for the turbulent quantities near the suction slot. 
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Figure 3: Contours of U velocity near the leading edge and over the hump. 
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